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MATHEMATICAL AND COMPUTATIONAL STUDIES OF EQUILIBRIUM 


CAPILLARY FREE SURFACES 

by Norman Albright, Nal-Fu Chen, 
Paul Concus, and Robert Finn 


We present here the results of our mathematical and computational 
studies of equilibrium capillary free surfaces that have been carried 
out over the past few years. These studies center on the equation 


( 1 ) 



= Ku + 2H , 


, ■ . 1/2 

W = [1 + (Su/ax)'" .+ (3u/ay)^ ] 


which is satisfied by the single— valued height u(x,y) of a liquid-vapor 
free surface in Laplace equilibrium (mechanical balance between 
gravitational and surface forces) . The gravitational force is assumed 
to be uniform and directed vertically, and in (1) the- quantity K = pg/cJ 
is the capillary constant, with p the difference in densities between 
the vapor and liquid phases, g the acceleration field — positive, 
negative, or zero — and cr the liquid-vapor surface tension; 2H is a 
constant, which is determined for a given problem, usually indirectly, 
by constraints such as container shape, liquid volume, and contact 
angle . 

Our studies have pursued several directions, each of which is 
discussed independently in a self-contained manner in Appendices I to V. 

In Appendix I, the general question is considered of whether a 
wetting liquid always rises hi^er in a small capillary tube than in 
a larger one, when both are dipped vertically into an infinite 
reservoir. For this situation, the height of the free surface 
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satisfies (1) , with H = 0 when u(x,y) is measured vertically upward from 
the reservoir level u=0, and the contact angle condition 

(2) I -S = 

is imposed at the tube walls, where 8u/8n is the derivative of u with 

respect to the outer directed normal to the wall and y is the contact 

angle (between 0 and 7 t/ 2 for a wetting liquid). 

It is generally known that if a given circular tube T^, dipped 

into a reservoir, is replaced by a concentric one of smaller diameter, 

then the liqviid-vapor free surface u^(x,y) in y is at all points 

higher than was the (original) surface u^(k,y) in T^. We have proved 

that for circular T this is still the case, no matter what the place- 

o 

ment or cross section of T^, provided only that it lies interior to 

the space originally occupied by T^. 

When the larger tube T is not circular, the situation can be very 
■ o 

different. In chapter I it is proved, by example, that if there 

exist configurations for which the volume (and hence the average rise 
height) of liqviid lifted by over the cross section of y exceeds 
that lifted by For the configurations considered in Appendix I the 

cross section of the larger tube -has a corner, however the result can 
be obtained even if comers are rounded so that both tubes are smooth 
and convex. We are continuing to pursue this question beyond the 
initial Investigation reported here. 
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In Appendix II is initiated an analytical investigation of the 
relationship between the family of solutions of (1) corresponding to 
rotationally S 3 nnmetric pendent liquid drops and the singular solution, 
corresponding to an infinite spike of liquid extending downward to 
infinity, whose existence we have proved previously. For these 
surfaces, K < 0, and H = 0 for u measured from the free plane reference 
level, so that the rotationally symmetric form of (1) , suitably 
normalized, is 


1 

r dr 


r du/dr 

r 91^/2 

1 + (du/ dr) 


-u 


9 


or, in parametric form, 


(4) 


iL 

ds 

du 

ds 

ds 


cos ij; 

sin ijj 
sin 

-u ^ 


Here s is arc length, x is the distance from the axis of sjrmmetry, and 
if) is the angle made by a tangent to a vertical section with the r axis. 
For the pendent-drop family of solutions, the initial conditions 

are 


u(0) = < 0, 


9 


(5) 


r(0) = 0, 


i|) (_0) = 0 
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whereas for the singular solution u ~ -l/r as r ^ 0. Fe prove in 
‘Appendix II that for sufficiently siaall in magnitude the solutions 
of (4,3) can be expressed as a sin^e-valued function u. = u(r) for all r, 
i.e. no verticals can occur. ' ‘ 

For larger values of ju^| we prove that no such global representation 
is possible. (Numerical computations give u^ “ -2.57 as the critical 
value for which a vertical first appears.) 

For very large values of |u^| the solutions of (4,5) are shown 
to have near their vertex the general fom of a succession of circular 
arcs joined near the vertical axis by small arcs of large curvature. 

J 

However, their continuation eventually projects simply on u=0 and 
continues in an oscillatory manner to infinity. Our analytical estimates 
for large |u^| are obtained using a technique that compares the solution 
surface with surfaces of constant mean curvature that are generated 
by the roiilades of an ellipse, . . 

We conjecture that as |u^j “ the pendent-drop solutions converge, 
uniformly in any |u| < u < “*, to the signular solution, which has a 
simple projection on u = 0. If the pendent-drop solutions are considered 
from the point of view of their global embedding in the manifold of 
all formal solutions of the equations, then the vertex of the drops 
appear as a transition point mar'king a change in qiialitative 
appearance. It is conjectured that the only global solution of (4) 
without double points, in this extended sense, is the singtalar 


solution. 
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Our study reported dn Appendix II has been extende'd , and will appear 
in Its completed form as reference [18] of that chapter. These studies 
of equilibrium pendent liquid drops are intended as a first step in 
approaching the question of stability of such drops. 

The studies in Appendices III and IV are computational ones that 
center on developing and comparing numerical methods for solving (1,2) 
for k: > 0 in domains that are noncircular — i.e., finding equilibrium 
free surfaces in vertical cylinders with noncircular cross sections 
for the case of downward-acting gravity. (If K 0, existence or 
taniqueness may fail.) Without a loss .of generality, only wetting 
liquids (0 ^ Y fr/2) are considered. Of particular interest are the 
difficulties that can arise for the case of a domain containing corners. 

The model problem considered in Appendix III is that of a cylinder 
with a square cross section. Because of symmetry, only one-quarter 
of the square is taken as the domain for the study, and we place on 
the domain a uniform square mesh, obtaining on it a corresponding 
discretization of (1,2). We investigate extensions of the block 
successive overrelaxation-Kewton method and of the- generalized 
conjugate gradient method for obtaining the solution of the discrete 
problem. Both methods are found to be satisfactory for obtaining 
accurate solutions for K > 0 for all cases in which the solution of 
(1,2) is bounded, although, at least for the model problem, the generalized 
conjugate gradient method is found to be more efficient and robust. 

As we have proved earlier, if a+y <71/2, where 2a is the interior 
angle at a comer of the cylinder cross-section (for the square a = 77/4), 
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the solution is unbounded at the vertex, the asymptotic behavior near 
the vertex being given by 

9 9 1/7 

(5) v(p,6) = [cos0 - (k -sin 0) ]/(klcp) , k = sina secy , 

where (p,6) are polar coordinates centered at the vertex, with 0 
measured from the interior angle bisector. Thus special numerical 
methods are needed in this case to obtain an accurate numerical solution. 

The methods we investigate in Appendix III are based on deleting a 
neighborhood of the singular corner from the numerical domain and using 
(5) in this neighborhood. The normal derivative of the asymptotic 
solution v(p,9) provides the required boundary condition for the 
numerical problem, and the resulting solution height is then matched 
at the interface. In choosing the position of the Interface one must 
strike a balance between keeping it close enough to the vertex so 
that the asymptotic representation is accurate over the deleted domain 
and keeping it far enough from the vertex so that the discretization 
in the numerical domain can represent the steep gradients near the 
corner adequately. 

In an attempt to overcome this difficulty, we investigate the 
possibility of overlapping the as 3 mptotic and numerical solutions over 
part of the domain. We find that such an attempt to improve the 
asymptotic solution by numerical means is not fruitful, in essence 
because of the steep gradients of the solution. 

From our experiments we conclude that the technique of deleting a 
small neighborhood of the corner from the ntimerical domain can succeed 
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in yielding a sufficiently accurate solution, particularly away from 
the comer. In fact, the solution over the bulk, of the interior of 
the square is rather insensitive to changes in the position of the 
interface between the deleted and mimerical domains. 

For reasons given in Appendix III,, we find it desirable (a) that 
the interface should be chosen to be a level curve of (3) (a circular 
arc intersecting the edges of the square with angle y), (b) that it be 
taken close to the singularity corner, (c) that a higher-order discretiza- 
tion, a mesh refinement, or other numerical method be used in the 
numerical domain, at least locally near the interface, to represent 
the steep gradients adequately. The conputer programs prepared in 
connection with the study in Appendix III were written primarily for a 
uniform square mesh and were not sufficiently general to permit 
inclusion of the above features . These are examined to some extent 
in Appendix IV. 

In Appendix IV, discretizations of (1,2) are solved on irregulary 
shaped domains using JASON, a general-purpose program for solving 
nonlinear elliptic equations on a nonuniform quadrilateral mesh. 

JASON is reliable and easy to use, but is not necessarily very efficient 
for .obtaining capillary surfaces, as it is designed for a very general 
• class of problems. It did provide means, however, for calculating 
capillary surfaces in more general domains than is possible with the 
experimental programs of Appendix III. With it, capillary surfaces 
are calculated on the ellipse, on a circle with reentrant notches.. 


and on other domains. 
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For the ellipse, comparisons are made between the calculated 
results and analytical estimates of the conditions for nonexistence of 
solutions of (1,2) on sufficiently curved domains in zero gravity. The 
calculated results are consistent with the analytical estimates, and 
indicate that the estimates are reasonably sharp for this case. 

For the circular cylinder with reentrant notches, and for a 
square domain, the case for which the solution surface is uriboiinded 
at a comer is studied by means of the technique described in 
Appendix III. Here, it is possible to approximate more closely a level 
curve of (5) for the interface between the deleted and the numerical 
domains. The numerical results for these cases are depicted graphically. 

In Appendix V, the analytical estimates for nonexistence of solutions 
of (1,2) on the ellipse in zero gravity are evaluated. It is found 
that if the ratio of the minor to major semiaxes of an ellipse is 
less than 0.6116, then a nontrivial critical contact angle exists such 
that for contact angles less than this critical one there is no solution. 
The critical angle estimates obtained here are the ones used in Appendix 
IV for the comparisons with the values observed from solving (1,2) 
numerically . 
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1 . In the “capillary problem” one seeks to determine a surface u(x) satisfying the 
equation 

div Tu = ku, Tu^J^Fu, W=T/l+~jFuj^, (1) 

in a domain Q, and the condition 


v-Tu — cosy (2) 

PS 

on Z = dQ. Here y is prescribed, k— is the “capillarity constant”, 

<7 

and V is the exterior directed normal on Z. The symbols p, g, a denote density, 
gravitational acceleration and surface tension (all assumed constant). For back- 
ground information see, e.g., [1,2]. In this note we assume the “contact angle” 
y is constant ^ and that ?c>0, corresponding to the familiar physical situation of 
liquid rising in a vertical capillary tube of homogeneous composition. 

M. Miranda has posed to us, informally, the question: let Qx Q, and let , w be 
solutions of (1, 2) in, respectively, and Is > u in ? 


2. The question is physically suggestive, in the sense one expects the rise 
height of liquid in a capillary tube of small section to exceed that in a tube of large 
section. 

If y=0 the question has an affirmative answer, under very general hypotheses, 
both on the domains Q, Qi and on assuniption of boundary data; this follows from 
the particular forms of the maximum principle given by Gerhardt [3] and by these 
authors [4]. 

A heuristic reasoning, using the methods of [4], will convince the reader that 
the answer is again positive in the case and I are spheres (not necessarily 
concentric). 

If y = Tcjl then the unique solution of (1, 2) in any Q is u(x) = 0; thus an affirm'a- 
tive answer (in an extended sense) is obtained in this limiting case, 

3 . In the present note we show that in general the answer is negative, eve^^^y^,^ 
convex domains with smooth boundaries. 

Consider the two dimensional region Q indicated in Fig. C. In § 7 of [5] is 
proved the existence of a unique solution u{x) of (1, 2) in Q. We .note the data y 


I 


We may assume 0 ^ y ^ f, as the complementary case reduces to this one by a formal transformation. 
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P. Concus and R. Finn 




need not be prescribed at the corner, where the boundary angle is not defined. 
At all other boundary points, the data are achieved strictly and smoothly [6, 7, 8], 
at least if y 0. 

We choose y so that a + y<f. Introducing polar coordinates centered at V, 
we set 


cos 0 — i//c-~sin^ 0 . sin a 

v{r,0) = ^ , k = . 

K k r cos y 


(3) 


There hold [4], 

div Tv = iCL: + 0(r^) 


(4) 


in Q, and 

V • Ta = cosy + 0(/-'^) 


(5) 


on the two straight segments. Further, there exists a constant C such that 

\u~v\<C (6) 


uniformly in a neighborhood of V. 

Let be as indicated in Fig. F. We note the circular arc T is a level curve for 
v; since ||7i-|>|i;J =: — v becomes infinite as r^O, there follows 

v-Tu],^! ' (7) 


as This relation is uniform on T. 

4. We apply the divergence theorem to i; in , obtaining 

IZj cosy— j V • Ti-'dcr^K J vdx + 0(\Z\^) (8) 

r fivOi 

by (4, 5). No contribution appears from the singularity at V, since |7y|< 1 for any 
' function /. 

Repeating the procedure with v replaced by u and taking the difference of the 
two expressions yields 

\\v ■ Tudo— ^ V • Tv da\^K J |« — u| dx + 0{iT|^). 

r Q I 


( 9 ) 
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Hence, by (6), 

J V- ri/r/(7>!r[ cosy - (10) 

r 

if r is chosen sufficiently close to V. 

5. Let »j(x) be the solution of (1, 2) in Oj . The divergence theorem now yields 

J V- Tiida-m cosy = ;c | (n-ujr/jc (11) 

r iii 

and we conclude, by (iO), that there exists an open subset of £2^ in which u>Uj . 
This completes the proof of our assertion, for the particular domains considered. 

6. An objection may possibly be raised, that the boundaries appearing in the 
above construction have singular points. However, an examination of our proce- 
dure, and of the method of construction of the function n(.x) in [5], shows that all 
corners — including the one at K— can be smoothed without affecting the qualitative 
result. The construction can also be effected so that £3jcc£2, and it can be re- 
peated without essential change in any number of dimensions. 

This work was supported in part by the Energy Research and Development Administration, by 
the National Aeronautics and Space Administration, and by the National Science Foundation. 

Note Added iti Proof. From Theorem 6 of [4] and known monotonicity properties of the spherically 
symmetric solution, one obtains an affirmative answer to Miranda’s question for any whenever I 
is a sphere. 
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Abstract 

A characterization is given -for the most general form of a symmetric 
pendent liquid drop^ including those situations in which equilibrium may not 
be stable. It is shown that for any vertex height u^ the vertical section 
of the drop can be continued globally as an analytic curve, without limit 
sets or double points. For small |u^| it is proved the section projects 
simply on the axis u= 0; . for large |u [ the section is shown to have 
near the vertex the general form of a succession of circular arcs joined near 
the vertical axis by small arcs of large curvature. However, the continuation 
of the section eventually proj ects simply on u = 0 and continues in an 
oscillatory manner to infinity. 

It is conjectured that as [u^j -> <» the section converges, uniformly 
in any |u| < u < <» , to a solution UCrj with simple projection on u = 0 
and an isolated singularity at r=0. The (liquid drop) solutions are 
studied from the point of view of .their global embedding in the manifold of 
all formal solutions of the equations. From this point of view, the vertex 
of the drop appears as a transition point marking a change in qualitative 
appearance. It is conjectured that the only global solution without double 
points, in this extended sense, is the singular solution U(r). 
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Tne form of the outer surface of a symmetric liquid drop suspended 
from a circular aperture is determined by the condition that the mean 
curvature of the surface is proportional to the distance below a hori- 
zontal reference plane. For points near which the surface can be des- 
cribed by a function u(x) we obtain an equation 



for the height u(x) above the plane. 

Here h is a physical constant, h > 0 when the liquid lies above the 
surface, and X is a Lagrange parameter, to be determined by the con- 
straints. 

In a specific problem the determination of X may lead to technical 
difficulties. Formally, however, X can be transformed out of (1) by 
adding a constant to u. in the present paper we intend to characterize 
all symmetric solutions for the case X=0, A solution corresponding to 
given X can then be found in this famiiy by transforming back. 

We shall also introduce the (inessential and convenient) normalization 
K - We then obtain, in terms of polar radius r, the equation 



for a symmetric two dimensional surface u(r). 

/ 

Not all surfaces that appesir pi^sicalfy have a simple projection on 
a base plane, hence for a complete description the form (2) is overly 
restrictive. We obtain a more suitable (parametric) form of the problem 
if we introduce the arc length s along a vertical section of the surface 
interface, measured from the vertex (0, u^). We are led to the system 
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dY 

ds 


- u 


1 

r 


sin Y 


(3) 


ds 


sin Y 


ds 


cos Y 


where Y is the angle between a tangent to the section and the r-axis, 
measured counterclockwise from the positively directed axis to the 
tangent line. 

From the point of view of general theory, one would expect a solution 
of (3) to be determined, at least locally, by the initial data 

{}) r(0) = 0 ; Y(0) = 0 ; u(0) = u© ; 

however, the system (3) is singular at s = 0, and because of this the 
second condition in (4) is superfluous (cf the discussion in [4] ), , 

The question of local existence has been studied by Lohnstein [5], 
who established the convergence of a formal power series expansion. 
Alternative ty, one could adapt the Picard method, as used by Johnson 
and Perko [ 6] for the capillary problem, to the case studied here. 

One obtains locally, by these methods, a non -parametric solution u(r) 
of the equation (2), which we may write in the form 

(5) (r sinY = -ru , 

corresponding to the (single) initial condition 

(G) u(0) =Uq. 


The circumstance that only one initial datum is required yields an 
impoi- cant- simplification for the problem of characterizing all solutions. 
It suffices to describe the. one-parameter family determined by u©, and 
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it is tnis approacn we adopt in ttie present work. 

In general, tne solution u(r;uo) determined in tnis way cannot be 

continued indefinitely as solution of (5). We snail snow nowever tnat 

for any Uq, tne function uCrju^) can be continued as a parametric 

solution of (3) for all s , yielding a surface witnout limit sets or 
double points . 

We snail cnaracterize quantitatively tne asymptotic form of tne surface 
in tne case ju^l :p 0, and we snail characterize qualitativefy tne global 
structure of all such surfaces. 

Tne global behavior changes qualitatively when [Uq! increases 
beyond a critical value. If | u ^ | » 0. there is an initial range for s in 
which the surface looks like a succession of spheres centered on tne 
u-axis with radius c2/ [ uj . In all cases, the section can be expressed 
for large s in the form u(r) and has an oscillatory behavior as r — > «o. 

If Uq = 0 the unique solution of (2) is given by u = 0. We assume 
throughout this paper that u q < 0 ; the remaining case is obtained by a 
simple change of sign We are interested particularly in what happens 
when uo<sc 0. Tne resulting surfaces are then physically unstable under 
most conditions of everyday experience; however, the problem has an 
independent mathematical interest (one specific feature of which we in- 
dicate below) and probably also a physical interest for situations in wnicn 
gravity forces are small compared with those of surface tension. 

We have proved in [7] the existence of a particular singular solution 
of (2) that can be expressed in the form U(r) in 0<r<6, and such 
that U(r) i as r — y 0. In [ 8] we have presented numerical evidence 
suggesting that the symmetric solutions discussed above tend uniformly 
to U(r) in any fixed region u>A>-oo. A particular consequence of tne 
analysis in the present paper will be a proof of a preliminary form of 
that conjecture, namely we .shall snow in section VI that the solutions 
converge asymptotically into a neighborhood of UCrj, the size of which can 
be estimated a priori and is small relative to other Clocal) distances. 
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In section VII the structure of solutions is examined from another point of 
view and some numerical results are presented, suggesting the types of possible 
global behavior. These results lead in turn to another conjecture, namely that 
the singular solution UCrj is the only global solution [in an extended sense) 
that is free of double points. . 

We remark that we know of few other studies of the problem from 

(3) 

a general theoretical point of view . To our knowledge the first 
attempt to characterize the shape of a liquid drop appears in Bashforth 
and Adams [ 10] in which a numerical procedure is developed to cal- 
culate the sectional form up to the first vertical point Thomson [ ll] 
used a geometrical method and was able to obtain a figure correspon- 
ding, in our notation, to u^cs -7. Computational studies were greatly 
facilitated by development of high speed computers and related techniques, 
and particular cases have now been calculated with much larger | u q | , 
see, e. g. , Hida and Miura [l2] and Concus and Finn [8]. Such calcu- 
lations are suggestive and instructive, but they cannot provide the 
unifying insight of a general formal description. The present work is 
intended as an initial step toward that objective. 

In this work we study the formal solutions of the equations and ignore 
the question of physical stabilily of the surfaces. With regard to this 
related matter the reader may wish to consult recent contributions by 
Pitts [ 13, 14] and by Hida and Miura [12], where also further referen- 
ces can be found. 

The central difficulty in the general study of the solutions of (2) lies 
in the failure of the maximum principle. In the particular situation 
studied here, a residue of this principle remains, permitting us to com- 
pare the solutions with those of a simpler equation. This circumstance, 
in conjunction with elementary formal manipulation of the equation, pro- 
vides the central tool in our investigation. We proceed in a succession 
of steps, most of which are elementary and immediate; when taken to- 
gether, however, they yield the requisite characterization. 

We remark that the comparison technique has proved effective also 
in other (related) contexts, and has led in particular to new information 
on the behavior of solutions of (2) near isolated singular points, see, 
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e. g., [15], 


The latter author wishes to thank J. Serrin and J. Spruck for a number 
of stimulating conversations. 


I The case of small i u q | . 


We shall prove : 

(4) 

Theorem 1 : If, in the initial value problem (5. 6) there holds 

Uq^ -2, then the solution can be continued as a (nonparametric) solution 

of the equation 



for all r > 0 , It has an infinity of zeros. For any two successive extrema 

r. , r,. of u(r) there holds 1 u{r^) I < | u(r ) | . Asymptotically as 

Uo — ^ 0 the first zero r is the first zero of the Bessel function J (r), 
o ■ o - ■ .i... ..I o 

r 2. 405. 
o 


We study first the portion of the -trajectory preceding the first zero, 
and we note that- (2) is equivalent to (5) on any interval on which [u^j <«®. 

li : Let u(r) satisfi^ (5) ^ 0 < r < R and (6) ^ r = 0. Then ' ' 
sin Y ( 0 ) = 0 . 

Proof : Integrating (5) from e > 0 to r, we find 

r 

r sin Y - e sin Y{6) = - / pudp , 

£ 

hence, using (6), 


(7) 


sin Y 


u. 


1 


1 + u 


2T 

r 



r 

/ p udp 
0 


from which we conclude lim u (r) =0. Hence there exists 

r* 

r 0 
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u(r) - u 

u (0) = lim ° 

r — >0 r 


1 

lim - / u (t) dr = 0 . 
r^O ^ 0 r 


lii 3Let u(r) satisfy (5) ^0< r< R and (6) a^r = 0. g u(r) < 0 
in 0 < r < R, then sin Y > 0 in this interval. 

The proof is contained in (7). 


It follows in particular that u(r) —r~^ 
sinY-D = li m sinY(r) exists, and that 

li 

r — ^ R 

1 ^ 

0<sinYj^=- g pu(p)dp ^1. 


Ujj^ ^ 0 as r 


R, that 


We conclude also that g the solution curve does not cross the hyper- 
bola ru = - 1, then sinYj^< 1. The following assertion covers as well 
the case of solution curves crossing that hyperbola. 


liii : Under the hypotheses of lii, g in addition u^> -2, then 
0< sin Y < 1 in 0< R. 

Proof : Consider the relation 
(8) + ( sinY = - u , 

the left side of which splits the mean curvature of the rotation surface 

defined by u(r) into a sum of latitudinal (k ■) md meridional (x ) 

.S. m 

sectional curvatures. We note lU that u(r) is increasing in 
0 < r < R; thus 


(9) 


sin Y 
r 


- - *72 


p u(p) dp > - 


u(r) 

2 


in that interval. Integrating (8) with respect to u and noting that 
Csin Y)^ = -(cos Y3^ yields, using (9), 

COS>Tj^>l- - 


which contains the assertion. We infer now from the general existence 
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theorem, applied at r - R, that the solution curve either can be continued 
upwaLrd until it crosses the r-axis, or else it tends asymptotically to 
this axis with increasing r. We may however exclude the latter possi- 
bility. 

r ^ 

liv : If u(r)< 0 in a< r< R<«, then R< aexpj- A— 

la sm Y 
L a 



Proof: From (5) we find r sinY > a sin Y in r< R. By lii, 
’ a 

sinY > 0. Thus, 
a 

, a sin Y 

du , a 

•r— = tan Y > sm Y > 


and the result follows on an integration. 


We have thus established that if u > - 2 the solution curve is in its 

o 

initial trajectory monotonically increasing and can be continued until it 
crosses the r-axis at a point r = a^. To study the further trajectory, 
we observe that the curve can be continued at least locally across the 
axis as a solution of (5), and we compare its inclination at a given 
height h with the inclination of the initial branch at an equal negative 
height. 


Iv : If the curve can be continued monotonicalib^ to a height h above 
the r- axis, then its inclination at this height is smaller than the incli- 
nation of the initial branch at the height - h , that is. 



dr 


h * 


Proof: We integrate (8) with respect to u between the height -h 
and h, obtaining 


- cos ^ 

- h 



sin Y 


cos Y 

h 


r 


du > 0 . 
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Ivi ; Under tne conditions of liv, the curve is strictly convex 

downward when u> 0, and u < - u. 

rr 


Proof: From (8), 


( sin y ) = 
r 


u 


rr 




3/2 


= - u - 


sin Y 
r 


< - h , 


From Iv and Ivi we find 

Ivii : The curve can be continued to a maximum height h, < |u I 
at a point r = m^ > at which point sinyCm^) = 0. 

We now proceed as above, comparing inclinations at corresponding 
heights until the curve crosses the r-axis a second time, then com- 
paring inclinations as in Iv, and so on. We obtain the qualitative pic- 
ture indicated in Theorem 1, of a curve oscillating about the r-axis 
with successively decreasing extrema (see Fig 1 ). We note also the 
additional information, yielded by the method: 


I viii : All inflections of the curve occur on (monotone) curve seg - 
ments approaching the r-axis, in the sense of Increasing s. At any two 
successive points a, 3, at which |u 1 = tu_|, there holds 


du 

> 

du 

dr 

a 

dr 


To prove the final statement of Theorem 1, we note (7), lii and 
Iviii that | u^(r;u^) | tends uniformly to zero with u^; thus the function 
v(r;u^) = u’^ u(r;u^) tends uniformly to J^(r) as u^ — > 0. 


II Large [ u^ j ; initial arc 

If h^« 0 the above reasoning on the behavior of the initial segment 
fails, and so do the results. 
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(4) 

Ti^eo££m 2i If ^o~ ~ ^ * there exists a value , beyond which 

u(r) cannot be continued as a solution of {5}. ^ r — sinY— 

The .proof could proceed by a direct study of the equation, as in 
section I. We obtain more precise results and also develop techniques 
that will be needed later if we proceed instead via an obvious comparison 
principle. 

II i : I,et v^^^(r), v^^^(r) be fvmctions defined in a< r< b and such 

that (rsinY^^^) > (rsinY^^^) . Suppose sin Y ^^^(a) > sin Y^^^(a). 

I — r r — = — 

Then sinY^^^(b) > sinY^^^(b), and, equality holds if and only if 

( 1 ) ( 2 ) 

v = V + const, on as ri< b.. 

The interest in II i lies in the fact that ^ (r sin Y is exactly 
twice the mean curvature of the rotation surface defined by u{r), and 
this circumstance facilitates the choice of comparison surfaces. In the 
present case we choose as initial comparison surface the sphere of 
constant mean curvature - u^/2, with center at the point 
(r, u) = ( 0, - 2/u^). Thus, if v{r) describes a vertical section of 

the sphere, there holds u(0) = v(0), u(r)< v(r) in the interval 0< r< - 
( see Fig 2 ). Using I ii, we find: 

II ii : The solution u(r) of (5, 6) can be continued at least until 

u r 

r = - 2 /u , and sin Y ( r) < - —r — . 

O -r~~ Z 

We need also; 

II iii : A solution u(r) of (5) admits' no Inflections in the region 
r u< - 1. 

Proof; From (8) follows ru+ sinY = 0 at any inflection. 

Thus, Y must continue to increase until either a vertical point is 
reached or the curve meets again the hyperbola ru = - 1. Integrating 
(8) with respect to u axid using II ii yields 


CJ)m 
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1 - COS Y> - ^ (u^-uu)- 

2 O 


From this we conclude that a vertical slope appears at a value 

(10) u^< (1+ 1 - 8/u^ } 

which completes the proof of Theorem 2. 

We may use a similar procedure to estimate the value We note 
that if w(r) describes a vertical section of the sphere of constant mean 
curvature 

(11) . i . . ( 1 + ]i- 8/u= ) 

with center at (0, + p ), then there holds u(0) = w(0) , and by II i 

u' (r)> w'(r) on any interval 0< r< R along which pu^ - 2. This 

condition is however satisfied at u * u^ (10), hence on the entire 

arc u <u<u,. We conclude u(r) > w(r) until the first vertical occurs 
o 1 

at r = r^ < P . 

We note that at r = p, where w'(r) -oo, the circle w(r) intersects 
the hyperbola ru = - 2. 

From II iii we conclude the initial solution curve is convex in the 
region ru< - 1. This property holds in fact for the entire arc; on the 
segment of u(r) joining the initial point to the point (r^, u^) on the 
hyperbola r u = - 1 we obtain from (7, 8^ I ii), using the comparison 
circle v(r), 

(sin¥)^= - u - > -■'+ > - v(r^)+ -^ > 0. 

The last relation holds whenever u ^ which is the condition 

o ’ 

that v(r) and the hyperbola ru--l, u<0 intersect*- We conclude 
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also from I ii and the relation 


(12) 


sin Y = - 


ru 1 i. 2 , 

~ + 27 i P 


that ru> - 2 on the arc considered. 


It turns out the sectional curvatures h « and h are both monotone 

Jl m 


decreasing on the initial arc. We have 


(13) 


d d sinY 2 / ^ u 1 / „ 

= rs / pudu -- = -(. 2 


sin Y 


dr dr 


- n}< 0 


0 


by (8, 12). Also, we have from (7, 8, 1 ii ) 


(14) 


^ K = ( sin Y ) = - u + - ( 2 + u ) 

dr m rr r r r 


1 * '1 ^ 

<-u + — (u-u )~- — / p u dp < 0 


r r 


0 


PP 


by the convexiiy of u. 


At the initial point (0, u^) there holds 
(8, 12, I ii) follows for r> 0 on .the initial arc 

(15) 

From (7, 8.i we have also 

u 

f-, sinY , o 

' ‘ m r 2 


- u /2 . From 
o' 


The inecualiiy (15) implies > which is the meridional 

curvature cf the comparison surface v(r). Comparing the surface u(r) 
with v(r) at corresponding values of u and applying II ii now yields 

(17) =“o- I 

o o 

(see Fig. 2) . 

We summarize the above results: 
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Theorem 3; Under the condition of Theorem 2, the initial arc of 
the solution curve^ from is convex, with sectional 

curvatures h , k* decreasing and satisfying x < -u/2<x* 
m X ^ m X 

in r < r r, . There holds 
— o 1 


u 


8 


- - < r,< - ( 1 - - -2 ) 

u 1 2 I 


(18) 


XL 

u - — <u<u - — (1 
o u 1 o 2 
o 


1 - ^2 ) 
u 

o 


For 0< r< - 2/u^ the arc lies below the comparison circle v(r) and 

has smaller curvature, and for u < u < u, the arc lies above the- 
— — o — 1 

comparison circle w-fy) and has larger curvature (see Fig. 2 ) . 


Further remark ; The hypothesis of Theorem 2 could 

he sharpened by using the comparison surfaces v(r) and 'w{r) in (7) 
and iterating. A direct numerical integration of (7) yields [l6] 
u^ « - 2. 5678 as the value for which a vertical first appears. We 
find immediately : 


II iv : Let be the largest value of u^ for which a vertical point 
appears. If u = u the vertical occurs at the second intersection of 
the solution curve -with the hyperbola ru = - 1, and is an inflection 
point for the solution curve ( see Fig 3 ). 


If u^S - 5, the upper bounds in (18) can. be expressed more simply. 


yielding 


(19) 


- 2 ^2 5 ^ 

< r, < - -- - — ^3 

u 1 u u 

o o o 

2 2 5 

u - — <u, <u - — -“3 

o u 1 o u u 

o o o 


These bounds could also be improved by iteration, starting with the 
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comparison surfaces v(r) and w(r). We note for reference that the 
asymptotic series obtained in [l6] by formal perturbation expansion 
yields, for the normalization used here. 


. ± - -^3 + 0 (u‘^ ) 

u 3u^ ' o ' 


( 20 ) 


u = u 
1 o 


u 


fLL^ + 0(u‘®) 

3 u ^ ^ o 

o 


as u. 


Ill Very large | u^j 

If u - 2 -^2 then u(r) cannot be continued beyond r, as a solution 
of (2). The curve can, however, be continued as a solution of the para- 
metric system (3) as long as r remains different from zero. We study 
now the behavior of this family of solutions in terms of the parameter u^, 

asymptoticalfy as u^ ^ - <». We base the discussion principally on 

II i; to do so, we introduce as comparison functions the sections of rotation 
surfaces generated Iry the roulades of an ellipse. The following result 
is due to Delaunay [ 17] : 

Let an ellipse of major axis 2a and distance 2c between focal 
points, roll rigidly on an axis without slipping . Let Q be the curve 
swept out by one of the focal points. Then the surface generated by 
rotating C about the axis has constant mean curvature H ~ (2a) 

We note that Q is periodic with half -period t satis^ing 2a< t< na, 
and that each half -period can be represented in the interval 
a-e<r<a+cbya Single valued function v(r) for which the equation 

(21) i(rsinY)^= 1/a 

holds, and for which sin y = 1 at the two end points (see Fig 4). 
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We proceed step hy step : 


The procedure of II shows that an infinite slope first appears at 
(ri, Uf), with bounds on given by (18). The system (3) is non- 

singular at (r u j), hence the curve can be continued beyond this point 
as a solution of (3,4). Prom (16) we find at 


H 


> 

m 


+ 





so the curve turns back toward the u -axis, and can be described again 

(locally) as a solution of (5). We compare it with a roulade y^^^(r) 

^1 

whose mean curvature is - — and for which a. + c = r (see Fig 3). 

Z X 1 1 

Since v ^^^(r,) = -», II i yields u <v hence u(r)>v^^^(r) as 
r 1 r r 

long as the continuation of both u and as single valued functions 
is possible. 


The curve v (r) can be continued toward the u-axis only until 

2 

the point (a^^ - c^> u^ + t^), with ~ “ y - r^ > 0 ; at this point 

^ 2 

the slope is again infinite. It follows there is a value r., > - — - r. 

Z Uj 1 

beyond which this branch of the solution curve cannot be continued as a 
single valued function. 


From the geometrical interpretation t ^ as the half-circumference 
of an ellipse with major axis 2a^ - - 2/u^ and focal length c^ - r^ - a^, 

one finds that for large 1 , 



Liet us estimate from above. To do so, we compare u(r) with 
a roulade Vj(r), which is determined by the conditions 





1 


u. 



( 23 ) 
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(23) ^ 

A » 2 2 2 ^ j A 

T. - f ja - c cos d© 

^ 0 

A formal estimate shows such a roulade exists if u^< ~ 2-^ . 


The conditions (23) are chosen so that the roulade can be placed 
with its lower vertical point at (see Pig 5), and so that in 

that configuration its mean curvature will be exactly the one determined 
from the right side of (5) by the upper vertical. Applying II i we obtain 

u^>v^^^^(r), u(r)<v^^^(r) for all r< r^^ for which u(r)<u^+x^. 

This condition clearly holds for r near r^; since v^^^(r)< u^ 
we conclude it holds on the entire interval < r < r^, thus 

0 > V ^^^(r) >u (r) > V ^^^(r) > -« 
r r r 

on this interval, and hence the solution can be continued to the left of r^, 
at least imtil the value 


(24) 


“"i ■ u, +¥ 


3 - -^ 1=^2 


For large ju^[ we find 


(25) 


with 


(26) 


^ _ 2 

T- + a. 

1 u^ 1 


ln|u,| 


u. 


A 

a 


1= -f ^ 0<snb ) 


Thus 

(27) 

with 


2 °-2 

^2 ^ " u, " 3 ■ 

1 u^ 


= 4 + ) . 


(28) 
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We now proceed, essentially, as in the proof of Theorem 2. We 
note 

sinY > sin ^ = - --- + “ ( 1 + ) ; 

2 r 2 


thus from (24-28) we find for r< ^ , 

sin Y ‘3 3,,^, -5,i 

” 50 “l + 0(“l ^i“ll )■ 


We integrate (8) in u from u(P ); using that cos Y < 0 until a vertical 
is reached, and that 

cosY02)> ) 

we are led to a contradiction unless the curve becomes vertical before 

-3 

u has increased by a value "16^^ . That is, a vertical must appear 

at a value 

(29) u^ < Uj + “ 16 u~^ 

The solution curve then turns back from the axis at {r 23 VL^ and initiates 
a further branch. 


We summarize these results : 


Theorem 4 ; From (r^,u^) the solution curve continues backwards 
towards the u- axis until a second vertical is reached , at a point (r„,u„) 
with 


(30) 


2 2 
— - r, < r_ < - -?r 

Uj 1 2 u^ + t j 


- = ^2 




( 1 ) 


In the interval r„ < r< r. there holds u < v 
2 1 r r 

in the interval 3„< r< r. there holds u > v 
2 1 r r 


( 1 ). 


u > V 

^ ( 1 ) 

U< V 

r 


We note in particular that the horizontal distance of the second vertical 
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from the axis exceeds that of the first vertical from the hyperbola 
ru = - 2. 

Ill i : There is exactly one inflection between (r^jU^) and (r^, u^). 
Proof; Clearly, at least one inflection appears. Using (8), we find 

( ru + sin Y ) = [ • ^ ■■ - ^ ) sin Y < 0 

r Vcos Y r / 

on the arc. Hence there is at most one inflection. 

We indicate briefly one further step in the procedure. We construct 

( 2 ) 2 
a roulade v (r) passing through with major axis 2ag = - ~ 

tn\ ^ 

and a second roulade v' ^(r) with a property analogous to that intro- 
duced for v^^^(r). Then there holds v ^^^<u < v 

r ■ r r 

in the intervals for which the comparison makes sense, and (as before) 
still another point (rg,Ug} is found such that sinY(rg) = 1. The proce- 
dure can be continued as long as the values of |u(r)| remain sufficiently 
large to justify the indicated steps. 

We find easily: 

III ii : The successive horizontal distances of the vertical points , 
from the axis and from the hyperbola, increase monotonically. 

Ill iii : On each arc segment returning from the hyperbola to the 
axis there is exactly one inflection. The same statement holds on the 
remaining arc segments for sufficiently large ju ] . 

Theorem 5: In the initial region u< 0, the entire curve is bounded 
(strictly) between the u- axis and the hyperbola ru = - 2 (see Pig 6). 

In this region, the curve can be represented by a single valued function 
r = r(u), with | r' (u) |. < <». 

Proof ; Since Theorem 4 and III i apply to any returning arc, we 
conclude the curve cannot contact the u-axis. The relation (9) shows 
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the curve does not meet the hyperbola on the initial arc starting from 

(0, u^) . To show this properly for any forward arc, we integrate (5) 

on such an arc from a vertical point (r u_.), obtaining 

2J 2] 


r sin Y 


23 


P 2 

r*. u . - r u 
23 23 


+ ^ / p^u'dp 


^3 


> - 


23 


2 

r u 




2 2 

from which^ru> - 2 on this arc. The same inequality shows sin Y > 0 
on this arc; an analogous integration establishes the same property on 
a returning arc. 


IV Global behavior 

The discussion thus far shows that the solution curve can be continued 
upward without self-intersections until it crosses the r-axis. For by I iii 
an outward branch must either achieve a vertical or cross that axis, 
and the comparison method of II yields readily that a returning branch 
has the same property. There are no horizontal points, by Theorem 5. 

We show here that a returning branch cannot cross the r-axis. Pre- 
cisely : 

IV i ; Let r = a^ be the first point at which the solution curve meets 
the r- axis . Then 0 < {a^^) <«». 

Suppose u'(a.-)< 0, or equivalently, cos Y - < 0. The curve could 

then be continued backward into the negative u-plane a first vertical 

(r ,u ) (see Pig 7), at which, by Theorem 5, 
a o, 

(31) r u > - 2 , 

a a 
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We integrate C®} with respect to Uj from to 0, obtaxning 

sin 'i' , m , ^ ^ 

(32) / du = cosY- + -r ^ 

' ' ti r 12a 

a 

To evaluate the left side of (32) we integrate (5) in r between r 

and r : 
a 

a 

r - r sin Y = - / pudp < 
a r 

2 



< 


r u 
a 

2 


+ r 

a 


by (31). Thus 


(33) > 


u 

g 

2 


on the entire arc. Placing (33) into (32) yields cos Y^> 0, contra- 
dicting the assumption. 


Now observe from (8) that at the crossing point a^^, the meridional 
curvature is negative; thus, if cos Y^ = 0. there would again be a back- 
ward branch from into the negative u-plane, and we obtain a contra- 
diction as above. 

From IV i one sees immediately that the proof of Theorem 1 applies 
without change to the region r ^ a^^, in the sense that the solution curve 
continues from the point (aj^,0) as indicated in Fig. 1. We show now the 
curve does not intersect the initial branch in the region u < 0. 

IV ii- Let u , Uo be two success ive points on the solution curve such 

T QJ 

that = Xg, with an intervening vertical at CryjnY)* then 

sin Yg < sin Y^; if . then sin Y^ > sin Y^. 

Proof: Suppose r^^ < r^. From (5) we find 

‘ 

Xg sin Yg - r^ = - / pu dp 

""a 

sin Y^-r,^ = - / pu dp. 
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thus. since x = x„, 
a 3 

r^Csin 'Fg- sin - / p(u“-u'*’)dp < 0 , 

u~ and u"*^ denoting values of u on the lower and upper branches. 
The case ^ follows similarly. 

From IV ii follows a < b in Fig. 8, and thus < b. But 

h. < h,, any j > 1, by I viii, hence h. < b, all j > 1, and 
1 3 

thus intersections are excluded. 

Combining these results with Theorems i and 5 we obtain: 


Tiaeorem 6; The solution of the parametric system (3, 4) defined by 
the data can be continued indefinitely as a non- self -intersecting 
curve. It has the form indicated in Figs 1,9, 10, 11. 


V Majgimum dia meter 

We define the diameter of a (symmetric) liquid drop as the largest 

diameter of all circular sections u = u., at -which the bounding surface 

■) 

is vertical. 
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From Theorem 1, 5, 6 we see that each drop has a well defined dia- 
meter. It is less obvious that there is a universal upper bound for the 
diameters of all possible drops, independent of u^. 


Theorem 7 : Let 6 ~ 2.473 be the unique positive root of the 
equation 

(34) r^ - 3^/^r - 3^^^ = 0. 

Then 2 6 exceeds the diameter of any solution of (3^4). 


We base the proof on a lemma, which also has an independent interest. 


V i : Let u(r) represent a solution curve passing through (a, u ) 
a 

with - 1 ^ au <0, and such that 
a 

(35) asin?^> a/2 , 

Suppose u(r) < 0 m a< r < R. Then sin Y > 0 on this arc segment. 

If the curve meets th e hyperbola ru = - 1 in a point (c, u ) with 

1 /4 

a<c<R, then c< 3 ' , and sinY > 1/2. 

Proof ; We integrate (5) between a and r, obtaining 

(36) r siny^ -a sinY^ = j ( a^u^ - r^u(r) ) + |-/p^u'(p)dp 

a 

from which, if a = a, 

I* 

(37) rsinY^^ - ~ r^u(r) + j / p^u' (p ) dp . 

3 . 

For r sufficiently near a, there holds sinY > 0. Thus, if siny were 

to vanish at any points interior to a< r < R, there would be a minimum 

r = r > a at which this occurs. But (37) would then imply 
^ r 

1 , Y 2 

0 = r siny > -r f P u'(p)dp > 0, 

Y Y a 

a contradiction. Thus, siny >0 on a< r< c, and hence u^(p)> 0 
on this interval. Setting now r = c in (37) yields 
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(38) 


sin y > 
c 


1 

2 


cu(c) 


1 

2 ' 


Finally, we note that at r = c the inclination of the solution curve can- 
not ceed that of the hyperbola. Thus, sin Y < - , and 

|1 + c 

4 

c < 3 follows from (38). 


We proceed to prove Theorem 7. For any given u^, the maximum 
width is attained at a point ^2j+l^ ^2j+1^2j+l ^ 

j> 0 (see section III ). At the preceding point (r ., u .) there holds 

2j 

either r„, = 0 (if j = 0), or else sinY-. = 1. In either event, (35) 

2j 2j 

holds with a = r .. Also - 1< r .u .< 0, and thus the curve crosses. 
2j 2J 4J 

the hyperbola ru = - 1 at a point (c, u^), < c< ^2j+l' 

a = c, r = r . in (36) and applying V i yields, using II iii, 
i 


(39) 




The (single) positive solution of (34) exceeds any solution of (39). 

Since j is arbitrary, ' we conclude 2 6 exceeds the diameter of any drop. 


VI Convergence t o the sin gular so lution 

The solutions discussed in this paper are apparently related to a 
singular solution U(r) of (2), whose existence we have proved in [7]. 

The function U(r) is defined in a deleted neighborhood of r = 0, and there 
holds asymptotically U(r) ~ , as r — >0. We have conjectured 
that in any interval 0< a^ r< b<<», the solutions of (3,4) admit 
a single valued representation u(r^u^) and converge uniformly to U(r), 
as u^ CO. Figures 9, 10,11 show the results of calculations supporting 

the conjecture. A preliminary (weak) form of the conjecture is proved in 
the complete, published version of this report [18] . 
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VII ISOLATED CHARACTER OF GLOBAL SOLUTIONS 

There is a strong numerical -evidence to suggest that global solutions 

of (3) lying in the region between the envelope of the solutions r(u;u^) 

and the u axis^ and without limit sets or self- intersections, are rare 

in the manifold of all solutions. We know of no such solutions, apart from 
this paper 

those described in and the singular solution UCr) . In Fig. 12 are shown 
samples of the result of numerical integration of (3} through the initial 
point p determined by the intersection of r(u;-8) with U(r3 and for 
varying initial angles a, measured counterclockwise from the arc of the 
curve r(uj-8) emanating from p in the direction of increasing u. 


We note the curves r(u;u^) can apparently be extended below the 
level u = u^, if the isolated (singular) point of contact with the u-axis 
is admitted. The point appears to mark a transition in qualitative behavior; 
above it, the curve behaves like a Delaunay arc generated by an ellipse 
(section IK). Below that point, the curve has the general appearanoe of 
a Delaunay arc generated by a hyperbola, with the characteristic double 
points of those arcs. 

An analogous transition occurs on neighboring solutions without 
occurence of singular points on the axis. In any event, if such singular 
points are admitted, the corresponding (extended) curves r(u;u^) are 
embedded naturally in a solution set, all of which develop double points 
for sufficiently negative u, with (we conjecture) the single exception 
of the solution U(r). 
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developing the material of Section VII and for preparing Figure 12. 


This work was supported in part by the Energy Research and Development 
Administration, by National Aeronautics and Space Administration, and by 
the National Science Foundation. 



42 


Footnotes 

I 

(1) pi For background information on the derivation of (1) see, 
e.g. [1,2,3], 

(2) p3 The remaining case can be realized physically, e, g, , as the 
lower surface of a column of water in a glass capillary tube. 

(3) p4 We call attention however to a remarkable existence theorem 
due to Wente [ 9] . 

(4) p5, 9, 12 This improves the result announced in [8], 

(5) p5 A stronger Result of this type is given in [4], 
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Figure 1 

The case > - 2 ; inflections 
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Figure 2 

Initial comparison surfaces , 


< -2 
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Figure 3 


The case u = u 

o oc 
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Figure 4 


Delaunay surface 
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Figure 8. Proof of Theorem 6 
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Figure 10 

= -8; singular solution 
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- - 16 ; singular solution 
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NUMERICAL SOLUTION OF THE CAPILLARY FREE 
SURFACE EQUATION ON A SQUARE 

Nal-Fu Chen* and Paul Concus 
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Berkeley, California 94720 

November 1977 
ABSTRACT 

Numerical methods are discussed for solving the partial differential 
equation describing the equilibrium free surface of a liquid in a vertical 
cylindrical container whose cross section has comers. Both the cases for 
which the free surface is bounded and for which it climbs infinitely high 
at a comer are considered. For the model problem of a container with 
square cross section comparative results of numerical experiments are 
given for a nonlinear relaxation method, for a conjugate gradient method, 
and for several techniques for handling the corner singularity utilizing 
the knoxm asymptotic form of the solution. Representative solution 
surfaces are depicted graphically. 


*Present address: ' Mathematics Department, University of Southern 
California, Los Angeles, CA 90007. ^ ^ 
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1. Introduction 

In this paper, we discuss algorithms for solTjing numerically the 
nonlinear partial differential equation describing the equilibrium 
free surface of a liquid under surface and gravitational forces. We 
are interested particularly in the case of a liquid partly filling a 
vertical cylinder, the cross section of which has corners. We assume that 
the surface height u(x,y) is a single-valued smooth function of x and y, 
that there is sufficient liquid to cover the cylinder base entirely, and 
that the gravitational field is uniform and directed vertically down- 
ward. Then the function u(x,y) satisfies 


9 r 1 ^ -Y , 

3x ^ W 9x •' 3y 


jr 1^ 8u _ , 

^ W 9y J 


2H, 


W = [l + (9u/9x)^ + 


( 1 ) 


where B = P^^g/cf with the difference in densities between the gas and 
liquid phases, g the gravitational acceleration (positive when directed 
downward) , and O the gas-liquid interfacial surface tension. The quantity 
2H is a constant determined, in general, by the shape of the cylinder 
cross section, the volume of liquid, and the boundary condition between 
the liquid surface and the cylinder wall. The mean curvature of the 
surface is H at points where u=0. We consider only the case B > 0. The 
parameter B is the ( dimens ionldiss) Bond number, if (1) is considered to 
have been made nondimens ional with respect to a characteristic length of 
unity. 

Let the prescribed angle of contact between the liquid free-surface 
and the cylinder wall, measured interior to the liquid, be denoted by 
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the constant y. This condition can he written as 

i 

W ^ ^ cylinder wall, (2) 

where 9u/9n is the derivative of u with respect to 

the outward normal. We consider here only the case of wetting liqixids, 
0^<7r/2. (The complementary case of non-wetting liquids can be derived 
from it directly by means of a simple transformation.) 

The value of the contact angle plays a crucial role in determining 
the qualitative nature of the solution when the domain under investigation 
contains comers. Let the interior angle of a corner be 2a; it is proved 

in j6,7j that the solution at the vertex of the corner is bounded if and 

< 

only if a + Y > ir/2. For the case a + y < set K^sin a sec y and 

introduce polar coordinates p, 6 with origin at the vertex and 0 
measured from the interior angle bisector. Consider the function 

v(x,y;y) = j^cosS - (k^- sin^0)'^ J /(kBp) , (3) 

for -a ^ 0 ^ a, p' > 0.' It is shown in j5j that there exists a constant 
C such that the solution u(x,y) satisfies ju(x,y) - V(x,y)j < C for 
sufficently small p. ‘Because of the differing qualitative natures of 
the bounded and imbounded solutions, we treat these two cases separately. 

We take the 2x2 square as the domain for our model problem. Because 
of the S 3 mimetry, we consider only the- equivalent problem on a 1 x 1 square 
with boundary conditions as shown In-'Figure 1. We place a uniform square 
grid of width h^l/^f on the domain and discretize (1) in a manner similar 
to that used in [3j for the minimal surface equation (B=H=:0) , which is 
related to the general procedure in Jll] for linear elliptic equations. 



W 



Fig.. 1 

The domain for the model problem. 
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We write the integral form of (1), as obtained by use of Green’s theorem, 
over a sub-domain D of which yields 


[Bu-b2H) dA, 


(4) 


where n is the outward unit normal and 3D is the boundary of D. 

'Then we place on an auxiliary mesh (the dotted lines bisecting 
the original mesh lines in Fig. 2). We denote an h x h square bounded 

by the original grid lines as a cell, and we denote a rectangle bounded 

’ i 

by auxiliary grid lines, and possibly 3f2, as an auxiliary cell. By 

imposing (4); on each of th^ auxiliary cells and using the midpoint rule 

I ' ' 

‘ ' » 

for integration, we obtain* a set of (nonlinear) algebraic equations. The 

I 

equation corresponding to a typical interior point u(x,y)=u(ih, jh)=U^ . 
is given by | ' ’ 


. . -W“^ f 2U. .-U. J .-U.. ! . ) i 

^>2 = ±,J\ I 

2U. .-U. . A 

1 2U. .-U. . ,-U. .._Y 


77 _ 

i+l,j’ 


+ W 


,^1 
i+l,jtl\ 


2U. . -D. .-U, , .T 1 
x,3 1+1,3 x,3+lj 


+2h { BU. .+2H = 0, 

1,1 


( 5 ) 




XBL77IH2I92 


Fig. 2 

The domain with the original grid and the auxiliary grid. 
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where 


w - ri+ 1 + n? T + n? , ■r')l 

I.J L \ i.J I.d-1 


with 


6 = (U. ,-U. , .]/h and n. t =(u. ,-U. . 

y 1,3 1-1,3/ 1,3 \ 1,3 1,3-1/ 


Notice that the integration of (4) is performed over each auxiliary 
cell, and we use the approximation / / (Bu+2H)dA-(BU. ,+2H)'(area of D) . 

JJd 

Note also that in the discrete equations we take ¥ to be constant over 
each cell =(Wt r for the cell centered at (i-l/2, i-l/2il . For a 


= (Wt r 
\ 1,3 


3-1/2|. 


boundary segment of we have 


/ 


an 


( —1 
Vu'ii ' cosy • (length of 9n^) , if W 3U/3n=cosY 

1 


W 


, if Su/3n = 0. 


In the following sections , we shall discuss and compare different 
methods for solving the discrete equations, for the bounded and unbounded 
cases . 

2. The Bounded Case 

In this section, we investigate two methods for obtaining the 
solution of the discrete form of (1) , (2) for the bounded' case 
(a + y ^ TT / 2)., which for the square (a - ^4) is the case y ^ ir/4. 
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2.1. Block SOR-Newton 


Here, we use the nonlinear block successive over relaxation method 
investigated in [4] for solving the mi nima l surface equation and 
described in [10] as , the one-step block successive overrelaxation- 
Newton (BSOR-Newton) method. A block SOR iteration, corresponding to 
a row of mesh points, is performed using a single Newton iteration to 
solve approximately the resulting system of nonlinear equations. We 
thus “solve" for one row of mesh points (see figure 2) at a time by 
considering the function values in that row U. 0 ^ i < N, to be 

the only unknoxms and using the latest values for all other U.. . . 

- 

If we denote by U_, f_ the vectors U_ = (U . U. , ..., U„ 

J d , J 0,3, x,3, N,j'" 

f_ = (f . f. .), and by the matrix 

J 0,3, 1,3 N,3 ’ JJ 

J ^ =(9f. ./9U .) ,i,k=0,l, . . . ,N evaluated at the current value of TJ. 

•J3 1,3 k,3 1, 

we have the iteration 




^(£+ 1 )=^^ ^(£) _ 
<J J 


where u is the relaxation parameter-. Jj.j is symmetric j positive-definite, 
and tridiagonal, thus the solution of ~ obtained efficiently 

J w xJ 

using Gauss elimination without pivoting, or Cholesky decomposition. 
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By differentiating (5) , we obtain explicit expressions for the 
partial derivatives at a general Interior point; 



and, by synmetry, 

=rd(w"^}/d!vu|l = - -y (l+U^+U^J“^^^ evaluated for cell i,j 


(see [3] [4]). 
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After each iteration we adjust the value of H, using (4) integrated 
over The left hand side yields 2 cosy and the right hand side, after 
discretization, BU + 2H, where U is the area-weighted average of U. .. 

The adjustment of H, which enhances the rate of convergence, corresponds 
to an adjustment of the volvime of liquid. A final adjustment can be 
made, if desired, after the solution has been obtained, by moving the 
solution vertically to correspond to a particular liquid volume. 

The convergence rate of the BS0R.-Newton method depends critically 

on the choice of the relaxation parameter co. Table 1 compares the 

niHnber of iterations required to decrease the residual !l f II 2 from its 

_6 

initial value of approximately 0.5 to ^10 , with zero initial approxima^ 

tion, Y = 7t/3, B=100, and h=l/40. 

As obtained from our experiments, the estimated as 3 naptotically 
optimal value. of o) for the above case is 1.58. In general, it may be 
necessary to use a value of w smaller than the optimal one to ensure 
convergence for a given initial approximation (see Tables 1 and 2) . The 
number of iterations required for convergence can be reduced in these 
cases by adjusting (O towards the asymptotically optimal value as the 
iteration proceeds, once the approximation to U becomes sufficiently 
good to permit doing so. Each complete SOR sweep (for the 41 x 41 grid) 
takes approximately 0.055 sec, on the GDC 7600 computer (using a FORTRAN 
program with the FTN4 compiler, OPT = 2). 
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We have tried an initial approximation of 2 .ero for this problem, 
and the one of U equal to the portion of the lower hemisphere satisfying 
(1) and (2) for B =0 and Y ~ TT/A. For these two initial approximations, 
the number of iterations required for convergence was essentially the 
same, when convergence- occurred. Changing the contact angle did not 
affect the -convergence rate substantially either', except that cpnv^ergence 
was more rapid when y was close to 'Tt/2. 

Table 2 . Number of iterations for convergence for- B=100, h=l/4Q, and 


\ 
Y N 

1.4 


1.5 

1.6 

tt/4 

72 


56 

divergence 

it/ 3 

72 


56 

44 


The size of B substantially affects the number Of iterations required 

for convergence. The smaller B is, the larger the nimiber of iterations 

« 

required to reduce the relative error below the desired tolerance, and 
the stronger the dependence on to. Table 3 compares behavior for a few 
values of o) for the problem, y = tt/ 4, B=1 with zero initial approximation, 
and ilf^^^Il 2 < 10"®. 
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Table 3'. 


Numiber of iterations for eotivergence for B=l, Y “ 'n~/4, h=l/40. 


U 



0 . 


dJlJ 

1.5 

1.8* 

1.84* 

number 

>200 

191 

123 

of 




iterations 





* u was set ecjual to 1.5 initially to prevent divergence and was 
later increased progressively to the indicated value. 


The slower convergence for smaller B should be expected, since for the 
boundary conditions (2), or those in Fig. 1, the problem becomes 
singular when B=0. 

Judging from the experimental results, if either we have a priori 
knowledge of an optimal value for w or if a scheme for improving o) is 
incorporated into the program, this overrelaxation method is quite 
efficient for larger values of B. We shoxild add here that this 
method works also for the case B~0, for which a closed-form solution is 
known. We discuss further experimental results for the behavior of 
the BSOR-Newton method in Sec. 2.3, where a comparison is made with 
results for the method discussed in the following section. 

2.2. Hybrid conjugate gradient method 

In this section, the conjugate gradient method is combined with 
a fast Helmholtz solver to obtain iteratively the solution to the discrete 
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form of (1) and (2). This method proceeds as follows (see [8]): 

(i) Given an initial approximation is a vector of length 

2 f-1^ 

(N+l) ). Arbitrarily define p^ . For k = 0, 1, 2, ... 


(ii) Compute 




and solve M? - -f . 

(iii) Compute the search direction 


pCK) ^ 

(fCk+i) (fcfl)] 

where Bk = )— 7 ^ — ^ 7 r-r ^ 0, m, 2m, 3m,, 

(f(fc)^ ^(h) ) 


^0 ^m " ^2m "• 


(iv) Compute the new approximation 


j(k)^ (H) 

k 

(f(k) ^(k)\ 

"V - (^(k)_ j^Ck)) 


and J is the Jacobian matrix (Bf . ./3U ] evaluated at ‘ . In our 

experiments we choose the restart parameter to be m = 9, as in [8], and we 
continue t-he iteration until EPS, where EPS is the desired 

tolerance. 

We choose the matrix M in step (ii) to be one that approximates J in 
some sense and for which a fast-direct method can be used to solve 

-f^^\ Our choice is the discrete Helmholtz operator scaled by a 
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diagonal matrix (see [ 8 ]). Specifically, M = where D 

- 2 

is a diagonal matrix whose entries are (9f^ j^^^i )/t2 diag (-A^)], 

and 2Aj^ is twice the discrete five-point Laplace operator on a mesh 

of width h, obtained by setting W El in the discrete form of (1) , (2) . We 

2 

make the choice K = 2Bh , so that M is identical to the discrete form of 

( 1 ) , ( 2 ) when WEI. 

(k) (k) 

To solve Mz = -f ' - we follow steps a to c below: 

(a) Compute 

(b) Use a fast solver to find x> where 

(^-2A^ + 2Bh^)x = 

(c) Compute = D 

(k) 

The multiplication of Jp in the calculation of ctj^ in step (iv) of 

the algorithm is carried out utilizing fully the sparsity of J. Since 

J is block tridiagonal and each block is itself tridiagonal the 

2 

multiplication takes less than 9(N+1) operations. 

In the table below, we list the number of iterations to obtain 
f^^ 2 for different y and B with h = 1/40, EPS = 10 and H E 1. 

Table 4. Number of it erations required for convergence for h = 1/40 

and U^°^E 0 . 



100 

1 

0.1 

0.01 


13 

22 

28 

40 

lr /3 

8 

10 

10 

10 
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Each iteration takes approximately .139 sec,, which includes .073 sec, 
for the fast solver, except for the first iteration, which requires 0.183 
sec. including preprocessing. The prograiB CMA (with parameter K=2) [2] was 
used to obtain the fast solution of Helmhpltg’s equation in our experiments. 
Notice that the dependence on the value of y in Table 4 is much stronger 
than it is for the BSOR-Newton method. The smaller the angle, the more 
nonlinear the problepij and the more iterations required because M is less 
good an approximation to J and hence to f. Here, the value of B also 
influences the number of iterations required for convergence, more so in 
the case y ^ 'n'/4, for which W becomes very large near the corner (0,0). 

The singular case B=0 can also he handled by this method when used with 
an appropriate fast solver, 

2,3, Comparisons 

We summarize -thi data from some of our experiments in the following 
tables. In all cases the initial approximation = 0 is used and the 

number of iterations given are those required to obtain a residual 



It appearg, from the data in Tables 5 and 6, that the hybrid conjugate 
gradient method perforing consistently better in terms of computer time 
for the model problem than does the block Qverrelaxation-Newton method. 

The conjugate gradient method has the further advantage of not requiring 
the estimation of an acceleration parameter such as [0 (the dependence on 
the value of the restart parameter m is not as significant), and it is 



Table 5 . Number of Iterations (CPU seconds) for the BSQR-Newton Iteration 


T~‘ '• 

t 



Estimated 

0) 

0) 

EPS= 







best 

for first 

after 20 

-3 

10 

10 

lG-=. 

lO"^ 

h 

B 

Y 

s 

to 

• O 

iterations 

iterations 


100 

ir/3 

1,58 

1.6 

1.6 

21(1.13) 

29(1,56) 

37(1.99) 

44(2.36) 


100 

ir/4 

1.59 

1.5 

1.55 

25(1.34) 

34(1.83) 

43(2.31) 

51(2.74) 

. 1 
40 

1 

tt/3 

1.84 

1.5 

1.84 

79(4.-24) 

98(5.26) 

116(6.22) 

>120(6.44) 


1 

7T/4 

1.84 

1,5 

1.84 

77(4.13) 

90(4.83) 

107(5.74) 

>120(6.44) 


.01 

tt/3 

1.85 

1,5 

1.8 

116(6.22) 

>120(6.44) 




.01 

tt/4 

1.85 

1.5 

1.8 

118(6.33) 

>120(6.44) 




100 

ir/3 

1.52 

1.2 

1.5 

15 (.2i) 

20 (.28) 

24 (.34) 

27 (.38.) 


100 

'rr/4 

1.52 

1,2 

1.5 

15 (.21) 

22 (.31) 

24(.34) 

27 (.38) 

1 

20 

2 

tt/3 

1.71 

1.5 

1.7 

46 (. 64) 

56(.78) 

66 (.92) 

76(1.06) 


1 

tt/4 

1.71 

1.5 

1.7 

46(,64) 

55 (.77) 

65 (.91) 

75(1.05) 


,.01 

5T/3 

1.81 

1,5 

1.8 

48 (.67) 

58 (.81) 

66(.92) 

• 79(1.11) 


.01 

tt/4 

1.82 

1.5 

1.8 

49 (.69) 

58 (.81) 

67 (.94) 

79(1.11) 









TABLE 6 . Number of iterations (CPU seconds) for the conjugate gradient iteration 


h 

B 

y 

EPS= 

10“^ 

10-"^ 

lO"^ 

10"^ 


100 

tt/3 ■ 

4(.60) 

5(.74) 

6(.88) 

8(1.16) 


100 

tt/4 

6(.79) 

8(1.16) 

10(1.44) 

13(1.90) 

1 

1 

tt/3 

5(.74) 

7(1.02) 

8(1.16) 

10(1.44) 

40 

1 

tt/4 

13(1.90) 

15(2.18) 

18(2.59) 

22(3.19) 


.01 

tt/3 

5(.74) 

6(.88) 

8(1.16) 

10(1.44) 


.01 

7r/4 

15(2.18) 

29(4.16) 

34(4.91) 

40(5.74) 


100 

tt/3 

4(.14) 

5(.17) 

6(,21) 

7(.24> 


100 

tt/4 

5(.17) 

6(.21) 

7(.24) 

9(.30) 

1 

1 

tt/3 

5(.17) 

7(.24) 

8(.27) 

9(.30) 

20 

1 

tt/4 

10(.34) 

14 (.48) 

16 (.54) 

19 (.64) 


.01 

tt/3 

5(.17) 

6(.21) 

8(.27) 

9(.30) 


.01 

tt/4 

13 (.44) 

22(.74) 

25 (.84) 

35(1.18) 
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less sensitive to the initial approximation, although it does require more 
computer storage than is required for the BSOR-Newton method [8] . For 
nonrectangular domains the conjugate gradient method would lose some of 
its competitive advantage of speed since more computer time is required 
by fast-direct methods in- this case for the solution of the Helmholtz 
equation. See [8] for other possible choices for M and for comparisons 
of the conjugate gradient method with the BSOR-Newton method for the 
case B = H = 0 with boundary conditions for which the problem is not 
singular . 

In order to estimate the discretization error, the numerical solution 
for B=0, Y = f/ 4, h = 1/20 was compared with the known closed-form 
solution. We found the relative error of the computed solution in the 
infinity norm to be less than about 1/2% everywhere beyond 3 grid 
points from the comer, 18% right at the corner, and 1 to 4% in between. 
The relative difference between the computed solutions on the grids for 
h = 1/40 and h = 1/20 was about 2% in the infinity norm for this case. 

3. The unbounded case 

As pointed out in Sec. 1, there is a critical contact angle y^, 
which is 7r/4 for our model problem, such that the solution is unbounded 
at the corner of the cylinder cross section for all Y 

asymptotic form of the solution in a neighborhood of the corner is given 
by (3). If Y < Yq> we use this asymptotic solution in conjunction with 
discrete methods away from the corner to solve our problem. 
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We investigate two methods that use different discretization proce- 
dures over the portion of the domain away from the corner. In the first 
method, a neighborhood of the corner is deleted entirely from the domain 
to be discretized.- In the second, the asymptotic solution is extended 
into part of this domain. We consider, as in the previous section, only 
the model problem of a square on which has been placed a uniform square 
mesh. 

3.1. Method A 

Here we assume that the asjrmptotic behavior (3) holds from the 
vertex of the corner with the singularity up to one or several grid 
points away. We term this portion of the domain the asymptotic domain 
and the remaining portion the numerical domain (see figure 3). 

A natural way to obtain the numerical domain is to cut along a 
level curve of the asymptotic solution. The analytical solution in this 
domain then can provide an upper bound on the solution u of (1) , (2) [9]. 
Although the level curves of the asymptotic solution are circular arcs, 
the limitations of our experimental computer program require us to 
approximate such an arc by a straight line passing through mesh points 
and making an angle of 45” with the edges of the square. Thus, we solve 
equations (1) and (2) numerically in the domain shown in figure 4. 

The boundary conditions of the previous section apply at all bound- 
ary points, except those on P, where we use a normal-derivative matching 
condition obtained by differentiating (3) . After solving the resulting 
problem in the numerical domain, we can obtain a value for v - u in the 
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domain 


Fig. 3 

Division of the domain' for Method A. 



Fig. 4 

The mimerical domain for the test 
problem. 


XBL77ll-2!88 
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asymptotic domain by matching at one (or more) points of T. In this 
way, we obtain an approximation to the solution over the entire domain 

We use the same discrete equation (5) at a general interior point 
of the numerical domain as in Sec. 2.1. On F, the discrete equation 
at a general point is (see figure 5) 
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and / di, is evaluated by Simpson's rule using the directional 

J -n W OlX 


T. 

1 


derivative from (3) for the values of W ^9v/3n, 


where 


W 8n 


a 

a 


1 

2 


-a^Cx+y) + a2(y-x) 

[2 (x^+y^) (l+aj^^+a^^)] 

- ^cos 0 -|k^ - sin^ (k B 

2 - 1 - 1/2 

f sin 9 cos 0 K - sin 0 -sin 0 


( 6 ) 


) 

y B 


( 0 , K, B, p as defined in (3)) . 

The discrete equation for points one mesh interval away from F 
in the interior of the numerical domain is the same as (5) , except 
that for (i, j) on r, _ is given by W_ = 






2„2 


. 1/2 


We consider only the BSOR-Newton method for solving the discretized 
equations in the numerical domain. General purpose programs for solving 
the discrete Helmholtz equation are being developed and were not 
available for use with the conjugate gradient method during our study. 
These programs require, in general, more computer time than ones for a 
rectangular^ ‘.domain. 

I 

The determination of how far into the domain F should be placed 
for optimal accuracy poses an essential difficulty for the method: If 
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r is too far into the domain, the boundary condition obtained from the 
as 3 nnptotic representation on T will be inaccurate; if F is too close to 
the corner, the discretization in the numerical domain near F'may not be 
able to represent well the steepness of the solution. Considering these 
two alternatives, it appears that it is better to put F close to the 
corner, because one can, in principle, always refine the mesh locally 
or use higher-order or other methods to handle the steepness of the 
solution, while errors committed by pushing F too far into the domain 
cannot be compensated for easily. 

In our numerical experiments, we choose F to be sufficiently 
close to the corner so that at p, the midpoint of F, W ^3v/9n > 

0.999 - cos 2.5°. We solve the resulting problem on the numerical' 
domain with a mesh size h of 1/40 and 1/80, ,to give an indication of the 
discretization error for the particialar choice of F. We also solve 
problems for nearby choices of F to determine the sensitivity to the 
positioning of F. 

Of course, the matching conditions (6) on F may not be accurate. 
Although the difference between u and the asymptotic solution V is 
bounded by a constant as the corner is approached, the difference between 
their derivatives may be large. If F is chosen to be an actual level 
curve of V, close enough to the vertex so that (Vv| is large and 
W ^8v/3n is essentially 1 there, then |Vu| would also be large and 
choosing W ^9u/9n to be essentially 1 along F should then be sufficiently 
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accurate to be in keeping with the discretization errors in the interior. 
However, along the straight-line approximation to the level curve required 
by our test program, W ^9v/3n can .vary appreciably, and attempting to 
match to W ^9‘u/3n as in (6) might lead to errors. Our goal is to 
obtain from our experimental program not necessarily solutions of the 
highest accuracy but rather an indication of the feasibility of our use 
of the asymptotic solution (3) near the comer singularity in conjunction 
with a discrete method elsewhere. 

We give here a summary of some of the typical behavior found in 
our numerical experiments. 

For the case y-Q° , B = l, and F at the position where 

W ^9v/3n[ =0.999 (i.e., F is on the 17th grid point of the first row in an 
P 

81 X 81 mesh, or the 9th grid point of the first row in a 41 x 41 mesh) , 

moving F one grid point in the 81 x 81 mesh produces less than .1% 

relative difference in the solution in the interior 90% of the domain, 
between .1% and 1% relative difference in a band closer to F covering 
about 7% of the domain, and between 1% and 3.5% for the points in the 
immediate neighborhood of F. The relative difference is correspondingly 
about twice as much in the 41 x 41 case when F is moved by one mesh 
point. The relative difference in the solution resulting from refining 
the mesh from 1/40 to 1/80, keeping F fixed, is less than 0.2% over the 

interior 50% of the domain, between 0.2% and 1% over 30% of the domain, 

between 1% and 5% over 10% of the domain, and between 5% and 15% on the 
points in the immediate neighborhood of F. 
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— 1 

For Y = 30®, -B = 1, and W 9v/8n[ =0.999, the relative differences 

P 

are about one-fourth of those of the above experiments for y = 0® . 

Less sensitivity is to be expected because the solution surface is 

generally not as steep in this case. For similar reasons, the case for 

Y = 0°, B = 10, W ■‘■3v/3n|^=0.999, has about only half of the relative 

differences found for the case Y “ B = 1. 

If we choose F so that W ^9v/8nj^=0. 9999, then the relative differences 

from refining the mesh were found tc be two to three times as large as 

they were for placing F at a position for which W ^3\)/9nj =0.999.> 

■“X 

Based on the above information, we choose F so that W 9v/3n| =0.999 

P 

for most of our experiments. 

For the, case of Y = B = 1, the solution height in the numerical 

domain was found to be Uj^=l,46 for the 81 81 grid. We compare 

solution heights obtained from placing F on the 9th grid point on the 

first row (in this case, W ^9v/3n| =0.9999) and placing F on the 17th grid 

P 

point on the first row (in this case, W ^9v/9n|^=0. 999) on a 81 x 81 
grid. Although rather large relative differences exist for points close 
to F, we observe that over 75% of the domain less than a 1% change occurs. 
This indicates that the solution over a large portion of the domain is 
rather insensitive to where F is placed. As mentioned previously, 
over the same large portion, refining the mesh from 1/40 ,to 1/80 also 
produces only minor changes. Thus, for that portion of the domain, 
a medium size mesh (say, h = 1/40) seems sufficient to produce acceptably 
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accurate results. We observe also that for fixed F the numerical solution 
appears to be converging as h tends to zero through the different mesh 
sizes. 

For the unbounded case, the convergence of the BSOR-Newton method 
is even more sensitive to changes of 0) than for the bounded case, 
especially for Y 0*^. For = 0, one must begin the iteration with 

(0 close to 1 to prevent divergence, and then quickly (within 20 iterations, 
say) increase w toward the estimated optimal value in order to obtain 
an acceptably rapid convergence rate. Overestimating w essentially 
always leads to divergence. For y = 0^ , the observed optimal values 
of CO found experimentally are listed in Table 7. 


Table 7 . Observed optimal co for y = 0°. 



The optimal values of (0 for Y ~ 30*^ about the same as for y “ 

One can see from Table 7 that both the grid size and value of B 
influence the value of the optimal (0. 

Table 8 summarizes our experimental results for a range of values 
of Y, B, and h. The number of iterations 



TABLE 8. Number of Iterations (CPU seconds) for method A 


Y B 


0 


0 

0 1 
0 10 
0 10 
0 10 




DEL = 


Initial 


after 


I 


approx. starting 20 steps (estimated) 


1.5 

1.79 

1.79 

1.5 

1.55 

1.55 


1.69 

1.70 

31(.46) 

42 (.63) 

51(.76) 

1.84 

1.84 

37(2.04) 

56(3.08) 

71(3.91) 

1.84** 

1.92 

70(16.1) 

112(25.76) 

170(39.10) 

1.59 • 

1.61 

24(.35) 

33(.49) 

41(.61) 

1.75 

1.78 

46(2.53) 

64(3.53) 

81(4.46) 

if 

1.79’^^ 

1.89 

90(20.70) 

141(32.43) 

178(40.94) 

1.69 

1.71 

24(.36) 

38(.57) 

48(.71) 

1.84 

1.84 

26(1.43) 

50(2.75) 

68(3.74) 

1.59 

1.61 

25(.37) 

34 (.51) 

42 (.63) 

1.75 

1.78 

48(2.64) 

66(3.64) 

84(4.63) 


— '1 ^ 

6 ^ 40 ^ 


10 ^ 1 


10 W 1 


The number of grid intervals from the corner at which F meets an edge of w. 
to = 1.9 after 100 iterations. 

0 - identically zero; S-portion of a sphere (i.e. exact solution for B=0, 'Y=it/4) 
to = 1.87 after 100 iterations. 


10 " 


61(.91) 
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given in Table 8 are those required for Ilb'^^^ll 2 ^ J5EL. 

A major difficulty for this method is the representation of the 
steep gradients of the solution in the ntimerical domain near F. In 
the next section, we investigate the possibility of handling this 
difficulty by overlapping a portion of the numerical with the 
asymptotic domain. 

3.2. Method B 

Here we investigate the possibility of Improving the previous 
method by overlapping the asymptotic and numerical solutions over part 
of the domain. ¥e divide the entire domain into the three regions 
shown in Figure 6. Region I is a small portion of the domain near the 
comer in which we assume that the asymptotic behavior (3) holds. 

Region III, including boundary 2 in Figure 6, is the purely numerical 
domain, in which we assume that the solution u is not too steep and can 
be computed accurately with the discretization used' previously. The 
region between is Region II, where we couple the asymptotic and numerical 
solutions under the assumption that u and v differ there by a lower- 
order quantity that is not steep and thus can be computed more accurately 
with a discretization than can either u or v* The staircase-shaped 
domain boundaries are used rather than the 45° lines of Method A- to 
approximate level curves because the limitations of our study permitted 
us to use only those boundaries that were simplest to include in the 
computer program. 

Let u be the solution to (1) and (2) , then in Region I and II we 



Numerical 


Boundary 2 


Boundary I 



Asymptotic 


XBL 7711-2189 


Fig. 6 

The three regions for Method B. 
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write u = V + e, where V is given by the asymptotic expression (3) and 
,|e[ is bounded by a constant. Over Region II, we consider e(x,y) = 
u(x,y)- V(x,y) and we derive discrete equations for e(x,y) involving 
the known quantity V, keeping only the first-order terms in e/v and in 
the derivatives of e divided by those of v. Over Region I, we take 


£ = e , as in Method A, where £ and its derivatives are negligible 

O O O D 

compared with V and its derivatives. 


In Region II, the left hand side of (4) becomes 


f 

[l + fv + + 


[ V + £ J 
y y' 


(v + e). 


^ 2 ^ 2 1/2 
+ + ^y 


-V £ -V £ 

2 2 
1 + Vy 




When we perform the integration along a direction parallel to y-axis, 
then the normal direction in (7) is parallel to the x-axis, and the 
right-hand side of (7 ) becomes (see figure 7) 



+ higher order terms. 
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If the normal direction n is parallel to the y-axis, as along C’, then 
the right-hand side of (6) becomes 


^ V dS, 


+ /£ 




2 2 
1+V/+Vy 


7 




+ higher order terms. 

To form discrete equations for e, we evaluate C7a) along C, for 

example, by using Simpson's rule xfith the five indicated node points for 

computing the three integrals involving V, with £„ replaced by 

(e.,, . - £. -)/h and e replaced by (e. . - - e, . + - e... ,)/21 

x+1,3 y 1,1-1 1,1 1+1, 3+1 1+1,1 

for the upper half of C and by (£■ . - e. . t + £.,, - ~ - i)/2h 

1,1 1 , 1—1 l "^!,! l **”!,!- 1 

for the lower half of C. If we then ignore terms that contain powers 
greater than one of e/V and of the ratio of their derivatives, we obtain 
a system of linear equations for ^ . The associated matrix is symmetric 
and banded. 


3.2.1. The it era t ion 


In solving the discrete equations in Region II, the boundary 
condition on boundary 1 is obtained from the directional derivative of the 
asymptotic formula (3), and the boundary condition at each iteration for 
boundary 2 is obtained by keeping the value of u fixed at its Region 
III value. We use the IMSL subroutine LEQTlB, designed to solve banded 
systems, to obtain the solution. Note that for the different boundary 
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conditions on boundary 2 arising from each iteration, only the right 
hand side of the system to be solved in Region II needs to be changed. 
Therefore, it is necessary to perform the associated triangular 

decomposition only for the first iteration. 

The method of solution for Region III is the BSOR-Newton method 
of the previous sections. The boundary condition on boundary 2 is that 
of keeping the value of U(i.e. e) fixed at its value obtained previously 
in Region II. Afterwards, we adjust the value of H in the same way as 
before, using the boundary of Region III as the contour for (3) . 

Each iteration thus consists of solving for e in Region II 
followed by solving for U in Region III and then adjusting H. 

3.2.2. Experimental Results 

We have experimented with this method for B=1 and B = 10 

on a 41 X 41 grid. We took boundary 1 to intersect the edge of the 
square one grid point from the corner, and boundary 2 to intersect the 
edge at the 10th grid point for B = 1 and the 4th grid point for B = 10. 

We found the solutions generally insensitive to these boundaries being 
moved by one or two mesh intervals. The time required for convergence 
was about three times greater than for Method A. 

Of importance is the observation, that in our experiments the 
value of |Ve[ obtained was not necessarily small compared with that of 
IVvJ. In fact, [Ve’l was almost half of [Vvj for most points in Region 
II, including the points close to the vertex. Thus the assumption, on 
which the numerical scheme in Region II is based, that |Ve| / [Vv| is 
small dges not hold eyen though. |ej /' |v[ ^ 0 as the vertex is approached. 
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3.3. Conclusion 

Since the results from Method B indicate that the assumption on 
which it is based does not generally hold, we conclude that Method B is 
not a suitable technique for improving Method A for this problem. The 
experimental results for Method A for the 21 x 21, 41 x 41, 81 x 81 grids, 
for a fixed position of P indicate convergence of the numerical solution 
as h(the mesh size) ->0. In addition, the values of U on the bulk of 
the domain change very little as the mesh is refined, indicating that 
reasonable accuracy can be obtained away from the corner, even for the 
coarser grids. ■ 

From the experimental data, it appears that taking F to be the 
line on which the value of at the midpoint is about 0.999 is 

reasonably satisfactory. Moving F one or two mesh Intervals produces 
less than 0.5% relative difference in the numerical solution over 90% 
of the grid points. As for the points close to F, one should consider 
using a finer mesh or higher-order discretization than in the rest of 
the domain and taking F to be a level curve of the asymptotic solution 
with an irregular mesh nearby. ^ These matters are considered further 
in a related study [1]. 

In Figures 9 to 18, we present graphical computer output depicting 
perspective views with contours for some of the numerical solutions 
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that were obtained with the methods of Section 2 for the bounded case 
and with Method A for the unbounded case. The perspective views 
displayed are those indicated in figure 8. The height contours are 

drawn at intervals of 0,1, measured from the center (1,1) .where the 

* 

height is taken to be zero. 

Solution surfaces for B=1 are depicted in figure 9 for 
(a bounded case) and in figures 10 to 13 for y= 30° and Y“0° (unbounded 
cases) . Those for B=10 are depicted in figure 14 for y= 60“ and 
figures 15 to 18 for y^SO** , All perspective views are at 

an inclination angle of 20“ and a viewpoint distance of 100 units. 
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viewing direction C^-) 




viewing direction (b). 


Figure 8 

Test problem domain showing perspective 
viewing directions. 





Figure 9 

Perspective view for B = l^from 
direction {a) , y = 60 . 



1 







3 
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Figure 14 

Perspective view for B = 10 from 
direction (a!) > Y = 60 . 





6 



1 
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Figure 18 

Perspective view for B = 10 from 
direction (b) , Y = 0** 
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ABSTRACT 

In this report we explore the nature of capillary surfaces on 
several noncircular ly symmetric domains and investigate the suitability 
of the computer program JASON for calculating these surfaces. We investi- 
gate how the accuracy of the computed solution depends on various parameters. 
We calculate capillary surfaces for elliptical cross sections . We describe 
a procedure for calculating surfaces on cross sections with corners for 
values of the contact angle less than the critical value. We use this 
method to calculate surfaces on a square cross section and on a circular 
cross section with a reentrant notch. 


■|/4 

intentiohalw Busriii 
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1. IimiODUCTION 

In this .report we explore. the nature of capillary surfaces on various 
noncircularly symmetry domains and investigate the suitability of the 
computer program JASON for calculating these surfaces. 

In the next section we define the capillary surface problem for 

liquids with contact angle y between 0® and 90®. The surface satisfies 

a pair of equations, Eqs. (1) and (2). In the following section we point 

out some elementary properties of these equations. Next we state some 

theorems first proved in Ref. 1' -for the nonexistence of solutions to 

Eqs. (1) and (2) in a zero gravitational field. If the curvature of the 

perimeter of the domain is sufficiently large at a point, then there exists 

a critical contact angle such that Eqs. (1) and (2) have no solution for 

Y < Y We give the value of y ... for a sharp corner and a lower 

‘ ' crxt ^ ‘ crxt 

bound on Y for a rounded corner. We then describe the nature of 

' crxt 

the solutions in a nonzero gravitational field. 

Next we describe JASON, a general purpose computer program that 
solves second-order elliptic partial differential equations in two dimensions 
on a nonuiiiform quadrilateral mesh by a finite element method. JASON is 
reliable and easy to use, but it is not necessarily the most efficient 
program for this problem, since it is designed for a very general class 
of problems. 

We consider the system of discrete equations corresponding to Eqs. 

(1) and (2) for a particular mesh. We study how the number of iterations 
and conputer time required to solve the discrete equations depend on the 
mesh spacing for uniform meshes. Next we explore how the accuracy of the 
computed solution varies over the cross section with the mesh spacing 



117 


and with the contact angle, on a square cross section for the case of 
zero gravity. A closed-form solution to Eqs. (1) and (2) is known for 
this case. 

Next we calculate capillary surfaces on elliptical cross sections. 

We investigate how the height at various points on the cross section 

depends on the contact angle and the ratio of the major and minor axes. 

Numerical estimates of Y . are obtained and compared with the lower 

crit 

bounds mathematically derived from a theorem given in Ref. 1. 

We describe a procedure for calculating capillary surfaces on cross 
sections with sharp corners for values of y less than We use 

this method to calculate surfaces on a square cross section and on a 
circular cross section with a reentrant notch. Finally, we discuss how 
the height of the surface at various points on these cross sections 
depends on the contact angle. 


2. DEFINITION OF THE PROBLEM 

We consider the equilibrium free surface of a liquid partly filling 
a vertical cylinder for the case in which the surface height u(x,y) is a 
single-valued smooth function of x and y and in which there is sufficient 
liquid to cover the base of the cylinder entirely. The gravitational 
field is taken to' be positive when directed vertically downward. The 
height then satisfies the equation 

V • Vu) = Ku + 2H 


W = (1 + |Vu|^)^ 


( 1 ) 
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where V is the two-dimensional operator (9/8x, 9/9y), K = pg/0 is 
the capillary constant, p is the difference in densities ‘between the 
liquid and gas phases, g is the acceleration due to gravity, and O is 
the gas-liquid surface tension. The constant 2H is determined by the 
cross-sectional shape of the cylinder, the volume of the liquid, and the 
boundary condition satisfied by the free .surface of the liquid at the 
cylinder wall. 

The boundary conditions for a free surface that makes a contact 
angle y with the cylinder wall is 

1 9u 

rr V" “ cos Y at the wall , (2) 

w 3n ' 

where 9u/3n denotes the derivative of u with respect to the outward 
directed normal at the wall. Only the case of wetting liquids, 0° ^ Y 90®, 
will be considered here. (The nonwetting case can be obtained from it 
directly by means of a simple transformation. ) 


3. ELEMENTARY PROPERTIES 

The quantity V * (W ^Vu) in Eq . (1) is twice the mean curvature of the 
surface z = u(x,y). Thus in a zero-gravity field (k = 0) the free 
surface has constant mean curvature. 

If one integrates Eq. (1) over the cross section of the cylinder, and 
uses Eq. (2), one obtains 

kAu + 2HA = L cos Y > (2) 

where u = A ^ J " u dx dy is the average height of - the free surface, and 
Si 

A and L are the area and perimeter, respectively, of the cross section 



119 


For the case K 0, if u(x,y) is a solution of Eqs. (1) and (2) 
for one value of H, then u(x,y) plus a constant is a solution for 
any other value of H. The solution is unique for a particular value of 
H if K ^ 0. The value of H is related to the reference level from 
which the height of the capillary surface is measured. It is sometimes 
convenient to measure the height of the surface from the level of an 
infinite reservoir into which the cylinder, without bottom, is dipped. 

On the surface of the reservoir, infinitely far from the cylinder, 
u = 0, Vu=0, and therefore H=0. The average height u^ to which the 
capillary surface inside the cylinder rises above the reservoir is then 
obtained from Eq. (3). 

u^ = (L cos y)/kA 

For the case (C = 0, one cannot arbitrarily choose H. It is 
determined by Eq. (3). 

2H =5 (L cos y)/A . (3) 

For IC=0, the solution of Eqs. (1) and (2) is determined only up to an 
additive constant. A convenient way to choose this constant is to take 
u = 0 at the lowest point of the capillary surface. For uniformity all 
the numerical results presented in this report will use this choice for 
u = 0, both for the case K=0 and for K 0. This will enable us to 
compare directly calculations for different values of K, 
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4. NONEXISTENCE OF SOLUTIONS FOR THE CASE K = 0 

For each cylinder, there may exist a critical contact angle such 
that Eqs. (1) and (2) have no solution for 0 ^ y < Several 

important theorems have been formulated in Refs. 1,4 that relate Y 

crit 

to the shape of the cross section of the cylinder. We give below forms 
of them relevant to our problem. 

Theorem 1 . Let the cross-section boundary have a corner with an 
Interior angle 2a. Then a solution of Eqs. (1) and (2) can exist in 
the neighborhood of that corner only if 

Y > 90“ - a 

and no solution exists if 

Y < 90° - a 

Example: a = 45° for a square cross section; thus the critical 

contact angle is y = 45°. There is no solution for 0 < y < y 

crit crit 

A similar result holds for a region on a curved boundary on which the 
curvature is large relative to a value determined by the entire cross 
section. 

Theorem 2 . Suppose there is a point p on the boundary at which the 
curvature is greater than L/A; then there exists a critical contact 
angle such that there is no solution of Eqs. (1) and (2) in a neighborhood 
oiB p for 0 < Y < Y _ . 

Example:- Let the cross section be an ellipse. The curvature of an. 
ellipse is greatest at the point where the semimajor axis cuts the boundary. 
The curvature at this point exceeds L/A. when the ratio b of semiminor 
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and semimajor ^xes is less than 0.6116. According ta Theorem 2, there 
exists a critical contact angle for an ellipse with b < 0.6116 [2]. 

A third theorem gives a lower bound for Y _ . Consider the cross 
section of a cylinder shoxjn in Fig. 1. Let A be the area of the entire 
cross section, and L be its perimeter. Let the cross section be cut by 
a curve T, which intersects the perimeter at points p^ and p^. The 
curve r divides the area of the cross section into two parts. A* and 
A -A*, and divides the perimeter into two parts, L* and L-L*. 

Let |r| denote the length of the curve T from to p^. Then 

the following theorem holds [1]. 


Theorem 3 . Let 


then 


V be the function 

V . ...JllA 

|l*a - a*l{' 

cos Y V 

crit o 


(4) 


where V is the minimum of V with respect to all possible curves F. 
o 

This minimization can be done in two steps. First, find the F that 
minimizes V for a fixed pair of intersection points, p^ and P 2 ; 
then vary p^ and p 2 * 

For any p^ and P 2 , the F that minimizes V is an arc of a circle 
with radius 


R 


LVp(p 




(5) ^ 


where V (p. ,p.,) is the minimum of V for fixed p.. and p„ U] . 
p 1 2 J- 

solution of Eqs. (4) and (5) gives Finally, minimization of 

V (p, ,Po) with respect to p^ and p„ gives V . 
p ± Z ± ^ o 
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Example: The lower bounds for ellipses with ratio b 

of semiminor and semimajor axes are: 


b 

Bound 

0.20 

35.12" 

0.25 

26.95" 

0.33 

16.88" 

0.40 

10.32“ 

0.50 

3.71" 

0.60 

0.12“ 


Numerical solution of Eqs. (1) and (2) indicates that these bounds are 


close to 


Y 


crit' 


5. BEHAVIOR OF SOLUTIONS FOR THE CASE K 0 

For the case K ^ 0, Eqs. (1) and (2) have a solution for each y. 

However, the behavior of the solution is different for Y Y 

’ crit 


Theorem 4 . Let the cross-section boundary have a corner with an 
interior angle 2a. Let p and 6 be polar coordinates centered at 
the corner with 6 measured from the bisector of the corner angle. Let 


- <^os Y _ cos y 
~ sin a cos 

where y = 90" -a, as in the case K=0. If y > y.., then u(p,6) 
crxt crxc 

is bounded in the neighborhood of the corner. 

If 0 < Y < then u(p,6) is asymptotic to 
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u(p,0) ~ u^(p,0) + [q cos 0 - (1 - q^sin^O)^] /Kp (6), 

as p-^0, where is 0(1). Inside the cylinder u grows like 1/p 

as the comer is approached. 

The solution in a cylinder with a curved boundary is bounded for 
all Y- However, this bound depends inversely on the radius of the 
corner [4] . 

Theorem 5 . Let C be the maximum curvature of the cross-section 
boundary, and let r = 1/C. Then in a neighborhood of the corner 

u(x,y) < r + ~ . (7) 


6. JASON 

The computer program JASON solves the nonlinear elliptic partial 
differential equation in two independent variables [5]': 

V • CcV(f)) = acj) + b (8) 

with boundary condition: 

cV4> • n = f(j) + g , (9) 

where a = a(x,y), b = b(x,y), c = c(x,y,|V^| ), f = f.(s), g = g(s), 
and s is the arc length along the boundary. It uses a bilinear finite 
element method on a nonuniform quadrilateral mesh. The resulting algebraic 
system of equations is solved by the one-step block successive over- 
relaxation-Newton method (BSOR-Newton) . 

The domain of a problem suitable for JASON can be composed of one or 
more distinct subregions separated by common boundaries . Each subregion 
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can have different functions a, b, and c. If the values of a, b, 
and c for a particular subregion are constant, they can be specified 
on a data card; if they are functions of x, y, or |V(|)j , JASON must 
be supplied with a subroutine to calculate them. One can suppress solution 
of Eq. (8) in a particular subregion by specifying c=0 in that subregion. 

JASON contains a mesh generator that produces a smooth mesh by 
"equipotential zoning." Input to the mesh generator consists of the set 
of boundary points of each subregion. Each point is specified by four 
numbers; the x and y coordinates of the point, and the (integer) L 
and M values of mesh line coordinates through the point. One need not 
specify every boundary point. JASON will compute by linear interpolation 
the boundary points that have been skipped, if they lie on the same mesh 
line. 

Dirichlet boundary conditions can be imposed at any point in the 

mesh. Neumann boundary conditions are restricted to the boundary of the 

entire domain of the problem and the boundaries between subregions with 
c = 0 and those with c 0. 

Before using JASON, one should make a rough sketch of the desired 
mesh, to guide the choice of the boundary points. These define a polygon 

whose boundary approximates that of the original domain. For the case 

K= 0, one must calculate the perimeter and area of this polygon, in 
order to calculate H. 

The input to JASON is a series of data cards. The first card 
specifies the procedures to be performed; generate a mesh, plot the 

I 

mesh, solve Eq. (8) for ^ , calculate V<f), and plot The second 

card is the title for printed and plotted output. The third card gives 
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the size of the mesh for the entire domain, LHAX and flMAX. The mesh 
coordinates take the values 1 < L < LMAX and 1 < M < MMAX. 

The next input is a set of cards for each subregion. The first 
card of each set lists the values of a, b, and c for a subregion, or 
provides parameters for a subroutine that calculates a, b, and c. 

The remaining cards of the set give the boundary points of that subregion. 

The next cards contain the values of L, M, f, and g at 
enough points to specify the boundary conditions. If f and g are 
constant along a boundary, one needs to give only L, M, f, and g at 
the first boundary point and L and M at the subsequent corners of 
that boundary. 

The next card allows one to choose the iteration control parameters 
of the solution subroutines . The remaining cards specify the domains 
for which the values of (j) or V<|) are to be printed or the domain in 
which $ is to be plotted. 

7. SQUARE CROSS SECTION WITH K=0 AND Y • 

The convergence properties of JASON were studied by calculating the 
capillary surfaces for ic=0 in a cylinder with a square cross section. 
The exact solution is known for 45° < Y 90^* ^nd is a portion of a 
lower hemisphere. The square had an edge of length 2 and was centered 
at the origin. It was the domain -1 x < 1 and -1 < y ^ 1. Figure 2 
shows the height of the capillary surface along, the edge of the square 
x = l for contact angles Y = 45°, 50°, 55°, and 60°. For Y = = 45° , 

the slope is infinite at the corner x = l and y = l, although the rise 


height is finite. 
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The same data are graphed in a different way in Fig. 3. This figure 
shows how the height at several fixed points on the edge of the square 
varies with the contact angle y. For y = 90" the surface is flat. 

Figure 3 shows that the height in the corner increases rapidly as Y 
nears The height at other points on the edge varies much less 

rapidly with Y. 

The coordinate axes x = 0 and y = 0 are lines of reflection 
symmetry, so the solution in one quadrant repeats in the other three 
quadrants. The solution was calculated in the first quadrant, 0 < x < 1 
and 0 y ^ 1. The boundary conditions are; n*Vu = 0 on the lines 
y=0 and x = 0, and (n*Vu)/W = cosy on the lines x = l and y = l. 

Here n denotes the outward directed unit normal for each edge. 

Capillary surfaces were calculated on uniform meshes with mesh 
spacing h = 0,2, 0.1, and 0.05. Surfaces were calculated for contact 
angle y = 45°, 50°, 55°, and 60° for each mesh spacing. 

Surfaces were also calculated using two nonuniform meshes, which 
were successive refinements of the uniform mesh with h=0.05. The 
refinements were made near the corner of the square. The first mesh 
(NUl) had h=0.05 for 0 < x < 0.8 or 0 < y < 0.8 and h=0.025 
for 0.8 < X 1.0 or 0.8 < y < 1.0. The second mesh (NU2) was a 
refinement of NUl which had h= 0.0125 for 0.9^x^l.0 or 
0.9 < y < 1.0. 

We now consider the number of relaxation cycles required to solve 
the system of discrete equations for a particular mesh. In subsequent 
paragraphs we will discuss how the solution of the discrete equations 
approaches the solution of Eqs. (1) and (2) as the mesh is refined. 
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The number of cycles required to solve the system of discrete 
equations depends on the mesh spacing. This was studied by calculating 
the surfaces for y = 50“ on a sequence of uniform meshes. Each mesh 
had n equally spaced horizontal mesh lines and n vertical mesh lines. 
The meshes were 6x6, 11x11, 16x16, and 21x21, JASON solves the 
discrete equations by the one-step block successive over-relaxation- 
Newton method, BSOR-Newton [6]. 

These results are compared with those of another program, which 
solves the discrete equations by a dynamic alternating direction implicit 
method, DADI [7] . Both computations were performed on a GDC 7600. The 
number of cycles and CPU-seconds required are shown in Table 1. 


TABLE 1. Comparison of BSOR-Newton with DADI 
on a square domain with n vertical and 
n horizontal mesh lines. 


Method 

Final 

R/L 


n = 6 

n = 11 

n = 16 

n=21 

BSOR-Newton 

1x10“^ 

Cycles used: ' 

22 

36 

63 

84 

BSOR-Newton 

ixio"^ 

CPU-sec used: 

0.615 

3.660 

13.152 

30.018 

DADI 

lxl0“® 

Cycles used: 

12 

16 

18 

24 

DADI 

ixio"® 

CPU-sec used: 

0.023 

0.-088 

0.204 

0.465 


The DADI method is approximately 18 times faster per cycle than JASON. 
However, JASON is a general-purpose program that solves the discrete 
equations^^ on a nonuniform quadrilateral mesh, and would be expected to 
be slower than the DADI program, which was written specifically for a 
uniform rectangular mesh for this problem. JASON is reliable and easy 
to use, but is not necessarily efficient for this problem. 
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The number of cycles required for JASON to solve the discrete 
equations depends on the criterion for terminating the iteration and 
the procedure followed in optimizing the relaxation parameter (i). The 
termination criterion used here is that R/L < 1 x 10 R is the 

square root of the sum of the squares of the residual at each point. 

L is the square root of the sura of the squares at u at each point. 

The initial value of 0 ) is 1.5. This was used until R/L was less 
_2 

than lxJ.0 , then (O was adjusted upward toward the estimated 
asymptotically optimal value, half the estimated increment every five 
SOR iterations. 

The number of cycles required depends on these iterations parameters . 
However, Table 1 gives a rough measure of the amount of computing necessary 
to solve the discrete equations, and how this depends on the mesh spacing. 

The number of cycles required to solve the discrete equations depends 
only weakly on the contact angle y* This can be seen by comparing the 
values of R/L for different Y after the same number of relaxation 
cycles. On a 21x21 uniform mesh, 80 cycles reduced R/L to 2.9x10 ^ 
for y= 45° and 2.1 xio for y~ 60°. On an 11x11 uniform mesh, 

30 cycles reduced R/L to 15.0x10 ^ for Y~^5° and 8.7x10 ^ for 
Y= 60“. This weak dependence on Y is shown in more detail in Table 2. 


^6 

Table 2. R/L x 10 after n cycles for various Y* 


Mesh 

Cycles 

Y = 45“ 

Y= 50° 

Y = 55° 

Y= 60° 

11x11 

30 

15.0 

12.5 

10.3 

8.7 

21x21 ' 

80 

2.9 

2.6 

2.3 

2.1 
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Table 3 shows an experimental study of the convergence of the solution 
of the discrete equations toward the known solution of Eqs. (1) and (2) 
as the mesh spacing decreases. The tabulated quantity 6u is the differ- 
ence between the known solution and the solution computed with a particular 
mesh. It is tabulated as a function of contact angle y, mesh spacing h, 
and position y along the edge of the square. The point y=l is the 
corner, and y=0 is the midpoint of the full square. 

Table 3 shows that for fixed values of y and h, fiu is larger 
for y near the corner. It further shows that for fixed values of y 
and h, 6u is larger for y near 45° . 

In the remaining paragraphs, we shall discuss how 6u changes with 
h for fixed values of y and y. We consider first the uniform meshes. 

Table 3 shows that for 0<y<0.9, 6u decreased by 1/4 when h was 

2 

decreased by 1/2- That is, ($u is proportional to h for these 

values of y on uniform meshes . However, 6u decreased less rapidly 
2 

than h at the corner y = l. Convergence is poorer at the corner 
probably because the slope is very large there. 

We consider next the nonuniform meshes. The mesh NUl has h=0.05 
for 0 < y < 0.8 and h= 0.025 for 0.8 < y < 1.0. We shall compare 
the 6u on MJl with that on the uniform mesh with h=0.05. For both 


the 

45° and 50° cases. 

6u 

decreased 

for 

0.9 < y 

<1.0, 6u was about 

the 

same for y=0.85. 

6u 

Increased 

for 

y =0.8, 

and 6u was about 

the 

same for 0 < y < 0. 

■6. 

That is, 

6u ( 

decreased 

in the upper half of 


the refined region. Increased at the boundary of the refined and unrefined 
regions, and was about the same in the interior of the unrefined region. 
The mesh NU2 "has h= 0.0125 for 0.9 < y < 1.0 and has h the 





TABLE 

3. fiuXlO^ 

for X = 

1.0 

and 

0 < y < 1.0. 




h 

y= 1.0 

0.95 

0,9 

0.85 


0.8 

0.6 

0.4 

0.2 

0.0 





Y = 

45° 






0.2 

34166 





2051 

487 

198 

108 

84 

0.1 

26689 


2114 



473 

102 

46 

26 

20 

0.05 

20944 

2239 

515 

161 


77 

25 

12 

7 

5 

NUl 

16761 

832 

220 

172 


140 

33 

11 

6 

2 

NU2 

13478 

434 

291 

217 


194 

29 

6 

3 

2 





Y = 

50° 






0.2 

4186 





831 

295 

143 

86 

70 

0.1 

1804 


503 



209 

73 

35 

21 

17 

0.05 

701 

■ 258 

125 

75 


51 

19 

9 

5 

5 

NUl 

304 

117 

81 

71 


64 

20 

8 

5 

4 


Y = 55° 


0.2 

1468 



391 

1^3 

94 


50 

0.1 

529 


188 

101 

44 

23 


12 

0.05 

176 

78 

48 

33 26 

11 

6 

4 

3 

Y = 60° 

0.2 

605 



192 

97 

57 

39 

33 

0.1 

199 


82 

50 

25 

14 

10 

8 

0.05 

62 

31 

21 

16 13 

6 

4 

3 

2 
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same as NUl for other values of y. "We shall compare the 6u on NU2 
with that on NUl. It decreased for 0.95 < y < 1.0, increased for 
0.8 < y < 0.9, and was about the same for 0 ^ y < 0.6. That is, 6u 
decreased in the upper half of the refined region, increased in the region 
of intermediate mesh spacing, and was about the same in the region with 
h= 0.05. 

8. ELLIPTICAL CROSS SECTION 

Capillary surfaces were calculated for cylinders of elliptical cross 
section with semimajor axis a = 1.0 and semiminor axes b= 0.60, 0.33, 
and 0.20. The cases b = 0.60 and 0.33 were calculated for K=0. 

The ease b = 0.20 was calculated for both K=0 and ic=l. Closed form 
solutions are not known for the ellipse as they are for the square for 
0 . 

Figure 4 shows how a quadrilateral mesh was fit to 1/4 of an ellipse. 

The 1/4 ellipse was made into a four "sided" figure by choosing a point 
on the curved edge and using it as a fourth corner. A mesh was fit to 
this by taking one pair of opposite sides to be mesh lines with constant 
mesh coordinate L, and the other pair to be mesh lines with constant 
coordinate M. 

Figures 5 and 6 show the 11 x 6 meshes generated for the cases 
b = 0.33 and 0.20. Figure 7 shows the 21 x 11 mesh generated for b = 0.20. 
Figure 8 shows the level curves of the capillary surface for the case 
b = 0.33 and y=60°. 

According to Theorem 2 there will be a critical contact angle if 
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the curvature at any point on the ellipse exceeds the ratio of its 
perimeter to its area. This occurs if b/a < 0.6116. Theorem 3 gives 
a lower bound for the critical angle. These are shown in Table 4. For 
the case K = 0, capillary surfaces exist for Y ^ ^crit’ 

Y < Y For K^O, .capillary surfaces exist and are bounded for y < 

but the bound is inversely related to both K and the maximum radius of 
curvature of the ellipse, as in Eq. (7). 

Theorems 2 and 3 restrict the existence of capillary surfaces for 
a cross section with a smoothly curved boundary. However, the actual 
cross section that is used is not an ellipse, but a polygonal approximation 
to an ellipse. Because of this. Theorem 1 gives an additional existence 
criterion, which depends on the interior angles at the polygon’s corners. 

By Theorem 1 the critical angle is “ 90° - a, where 2a is the 

smallest interior angle. The Y ... from Theorem 1 are shown in Table 4. 

Theorems 1, 2, and 3 limit the existence of solutions of the differ- 
ential Eqs. (1) and (2). However, the discrete system of equations 
approximating (1) and (2) can have a solution even when Eqs. (1) and (2) 
do not. This is worth remembering because the 11><6 meshes used to 
explore the ellipse problem are rather coarse, especially for the largest 
ellipse, b = 0.60. Thus the local truncation error, which is of the 
order of the mesh spacing squared, is not very small. 

The computing times used to calculate the various capillary surfaces 
are shown in Table 5. For the ic=0 cases, the smallest time was for ' 
the largest contact angle. For example, for the b = 0,60 case, with the 
11x6 mesh, 1.13 CPU seconds were used to calculate the surface for y = , 

and 1.57 CFU seconds were used for Y“0°. However, for the k= 1 case, 
the smallest time was for the smallest y. Note that it took approximately 
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TABLE 4. from Theorems 1 and 2 for various 

meshes and values of b. 


Y 

' crxt 


b 

Mesh 

Theorem 3 

Theorem 

0.60 

11 X 6 

> 0.12® 

8.79® 

0.60 

21X11 

> 0.12° 

4.55® 

0.33 

11 X 6 

>16.88® 

15.71® 

0.20 

11 X 6 

>35.12® 

24.89® 

0.20 

21X11 

>35.12® 

13.43° 


TABLE 5. CPU-sec used for various meshes and 
values of b . 


h 

K 

Mesh 

CPU 

- sec 

0.60 

0 

11 X 6 

1.13 

- 1.57 

0.60 

0 

21x11 

19 .,23 

- 19.39 

0.33 

0 

11 X 6 

3.66 

- 6.24 

0.20 

0- 

11 X 6 

6.20 

- 11.60 

0.20 

0 

21x 11 

44.70 

- 71.57 

0.20 

1 

11 X 6 

7.71 

- 9.38 
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the same time to compute surfaces for the ellipse with its irregular 
mesh as for the rectangle with its regular mesh. The rectangle with 64 
mesh points took 1.22 CPU. seconds. The ellipse with the 11x6 mesh 
took 1.13^1.57 CPU seconds for the b = 0.60 case. The b = 0.33 and 0.20 
cases took more time, but these surfaces had steep gradients. 

The computing time depends on the iteration control parameters. 
Because the 11x6 mesh is rather coarse, R/L < 1.0 x 10 ^ was used in 

most calculations for the criterion for terminating the iteration. For 

-5 

the 21x11 mesh a stronger criterion was used, R/L < 1.0 x 10 . The 

default value that JASON assumes is 1.0 x 10 

Three iteration parameters can be chosen to attempt the minimization 

of the computing time. They are; the initial value of (i), the value of 

R/L at which one begins to adjust u towards its estimated optimal 

value, and the fraction of the distance moved toward this value. The 

-3 

default values of these parameters are 1.5, 1.0 x IQ , - and 0.06. 

However, it was found that for most of the capillary surface calculations, 

JASON can solve the discrete .system of equations if the values 1.5, 

-2 

l.OxiO , and 0.60 are used. Computing time is less with these values 
because fewer steps are used before changing w and optimization is 
pursued more vigorously. 

JASON prints the values of R/L and other quantities at every tenth 
Iteration. If the solution 'attempt fails, one can examine this sequence 
of values and decide how to change the iteration control parameters. If 
R/L increases during the cycles when the initial w is being used, then 
the initial (o should be reduced, say to 1.2. If R/L decreases while 
the initial to is being used, but increases after optimization begins, 
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the fraction for the 0 ) optimizer should be reduced. If this does not 
succeed, it may be necessary to reduce the value of R/L at which 
optimization begins. 

Figure 9 shows the height of the capillary surfaces at the point A 
where the semimajor axis intersects the cylinder wall. Figure 10 shows 
the height at the point B where the semlminor axis intersects the 
cylinder wall. The height is taken to be zero at the center of the 
elliptical cross section. Figures 9 and 10 show u(A) and u(B) as 
a function of the contact angle Y for the four cases: b = 0.20, 0.33 
and 0.60 each with K=0, and b = 0.20 with K = l. Figure 9 shows that 
for fixed y, u(A) is larger for smaller values of the semiminor axis b, 
while Fig. 10 shows that for fixed y, u(B) is larger for larger values 
of b . 

We shall discuss first the case b = 0.20 with K=0. In this case 

the effect of the curvature of the actual ellipse near its major axis can 

be distinguished from the effect of the corner at point A of the polygonal 

approximation to the ellipse. The corner by itself would cause a 

of 24.89° for the 11x6 mesh and one of 13.43° for the 21x11 mesh. 

However, the ellipse's curvature overwhelms this effect and causes a 

y^rit ^ 35.12°. Figure 9 shows that u(A) becomes large with increasing 

rapidity as y approaches 35°. Furthermore, when the mesh is refined 

from 11x6 to 21x11 j the value of u(A) a^t 35° Increases by 12%, while 

that at 40° increases by less than 1%. This must be the effect of the 

ellipse's curvature, since the effect of the approximating polygon's corner 

is reduced for the refined mesh. This indicates that the y due to 

'crxt 

the ellipse's curvature is near 35°, and that the lower bound given by 
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Theorem 3 is indeed close to Y 

'crit 

Consider next the case b = 0.20 with K = l. The capillary surface 

for a smooth ellipse is bounded at A because it has finite curvature. 

However, u(A) should become infinite for Y less than the y , of 

crit 

the approximating polygon's corner. Figure 9 shows that u(A) for the 

11x6 mesh becomes larger as y approaches 25°, as is expected. Although 

the values of u(A) are quite different for the cases 1C=0 and tc = l, 

the values of u(B) are almost the same. At 35° they differ by only 0.3%. 

For the case b = 0.33 with K=0, the effects of the ellipse's 

curvature and the approximating polygon’s corner cannot be distinguished 

for the 11x6 mesh. The first predicts a y . ^ 16.88° and the second 

n ~ 15.71°. Figure 9 shows that u(A) and |du(A)/dy| become 

large as y approaches 15°, which is consistent with these predictions. 

It also shows that the due to the ellipse’s curvature cannot be 

much larger than 16.88°. This illustrates that the lower bound given by ' 

Theorem 3 is again close to y . . 

Finally, for the case b = 0.60 with tc = 0, the coarseness of the 

mesh caused the solution of the discrete system of equations to differ 

somewhat from those for the cases b = 0.20 and 0.33. The sharp corner 

should give a y of 8.79° for the 11x6 mesh and 4.55° for the 21x11 

crit 

mesh. The lower bound on v , due to the ellipse's curvature is 0.12°. 

'crit 

However, Fig. 9 shows that although u(A) increases as y approaches 0°, 
jdu(A)/dy| decreases. When the mesh was refined from 11x6 to 21x11, 
u(A) increased by 9%, 7%, 4%, and 2%, at y = 0°, 5°, 10°, and 15°, 
respectively. The relatively larger Increase for 0° and 5° than for 15° 
shows the effects of the corner, but not as dramatically as for the b = 0.20 
case. The coarseness of the mesh is shown by the change in u(B) when the 
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mesh was refined. It changed by 6%, 5%, 4%, and 3%, at y = 0®, 5°, 10", 
and 15°, respectively. In contrast, for the case b = 0.20 with K=0, 
u(B) changed by only 0.6% when the mesh was refined. 


9. CAICTJIATING CAPILLARY SURFACES ON CROSS SECTIONS WITH CORNERS 

A cylinder whose cross section has a corner, as in Fig. IIA, will 

have a critical contact angle. If the interior angle at the corner is 2a, 

then “ 90° -a. There are no solutions to Eqs. (1) and (2) for K=0 

and Y < Y - ^ . Solutions exist for K ^ 0 and y < J . ^ , but the height 
crit crit 

goes to infinity at the corner. It would be difficult to compute such a 

surface numerically over the entire cross section. Instead, the cross 

section is divided into two regions, as with curve AA’ of Fig. IIA. 

One is a small region Q that Includes the corner: the other 

c r 

is the remainder of the cross section. The asjnnptotic form of the surface 
near the corner region is knoxm. The surface on the remainder can be 

calculated numerically, with the boundary condition that it match the 
as 3 nnptotlc solution along the curve AA’ . 

The boundary data are the values of the function (n*Vu)/W at each 
point on the boundary. W is the function (1+ |Vu[ ) , and n is the 

outward directed normal at the boundary. (n*Vu)/W is cos y at the 
cylinder wall. It is zero on boundaries that are lines of reflection 
symmetry. It is matched with that of the asymptotic solution along the 
boundary curve AA' . 

The surface in a neighborhood of the corner has the asymptotic form [3] 

u(p,0) ~ u^(p,0) + v(p,6) 


( 10 ) 
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where is a function of 0(1) and 

v(p,0) = ^ [q cos Y - (1 - q^sin^ 0)^^^] . (11) 

The quantities p and 0 are polar coordinates of a point on the cross 
section* p is the radial distance from the corner, and 0 is the angle 
measured from the bisector of the corner angle. These are shown in Fig. 
IIB. The quantity q is 


q 


cos Y 
, sin a 


cos Y 


cos Y 


crit 


( 12 ) 


For p « 1, we have u « v. 

o 

The boundary condition function (n*Vu)/W can be broken into two 
factors, a magnitude and the inner product of two unit vectors. 

^ » « /S 

n* Vu _ I Vu I n*Vu , 

W “ W [Vu| ■ ^ 


We do not know the order of Vu - but we know that IVvl 

o ' ‘ 

so that in general jVu[ also must be at least 0(l/p^). 


MJ 

¥ 


= 1 - 


2lVu| 


= 1 - 0(p'*) 


is 0(l/p2), 
Consequently , 


(14) 


Thus if the curve AA' is sufficiently near the corner, then jVu[/W 
is approximately 1. 

We next consider the factor (n-Vu)/jVn|. The quantity Vu/[Vu| is 
a unit vector normal to a level curve of u. If we choose the curve AA' 
to be this level curve, then n = Vu/jVu|, and so (n«Vu)/|Vu| will be 1. 
The level curve of u is not known, but for small p we have u « v, 
so a level curve of u is approximately a level curve of v. We shall 
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choose the curve AA' to be a level curve of vj thus the factor 
(n*Vu)/|Vul will be approximately 1. 

The level curves of v are a family of circles, which depend on the 
contact angle y ^h® radius of a circle is 1/kv, and its center lies 
on the bisector of the corner at a distance q/Kv from the corner. The 
circle intersects the cross section of the cylinder wall at an angle y, 
as' shown in Fig. IIC. 

Let u denote the solution of Eqs. (1) and (2) on the entire cross 
section with (n*Vu)/W = cos y the cylinder wall. Let u^ denote 
the solution of (1) and (2) on with (n‘Vu)/W «= cos y at the cylinder 
wall and with (n*Vu)/W = 1 on the curve AA' . It is shown in Ref. 8 
that u^ is an upper bound to u at each point in 


10. SQUARE CROSS SECTION WITH ic=l 

The capillary surface over a square cross section was calculated for 
y = 30‘* and k= 1 by deleting the corner as described in the previous 
section. The value of ® square is 45®. The calculation was 

performed on 1/4 of a 2x2 square. The 9x9 mesh that was used is shown 
in Fig. 12. The boundary condition function (h*Vu)/W is 0 along the 
reflection symmetry boundary BOB'. It is cos 30° along the edges BA 
and B'A', and is approximately 1 on the curve AA' that divides the 
corner from the remainder of the cross section. Curve AA' is composed 
of eight straight line segments, which approximate a level curve of v 
for y = 30°. 

The heights of the capillary surface along the edge AB and along 
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the diagonal CO are shown in Fig. 13. The height is taken to be zero 
at the center of the square, 0. The height is approximately 1.80 on the 
boundary curve AA* and is almost constant along the curve. It differs 
by only 0.0062 between the endpoint A and the midpoint C of the curve. 

The point B is the midpoint of the edge of the full square. The 
height of the capillary surface u(B) is shown in Fig. 14 for K=1 and 

Y = 30°, 45“, 50°, 55°, and 60°. The four surfaces with Y ^ Y • were 

calculated on a 6x6 mesh with uniform spacing. In the case K=0, 
discussed previously, it was 'found that this mesh gives an accurate value 
for u(B), although the height at the corner is inaccurate for Y“'^5°. 

It is interesting to note that u(B) changes almost uniformly over the 
range 30° < y 90°. It does not change rapidly near Y=45“ or in the 

range 30° < y 45°. Although the fluid is drawn up to an infinite 

height in the corner, our calculation indicates that this does not affect 
dramatically the height at the midpoint. 


11. CIRCULAR CROSS SECTION WITH A REENTRANT NOTCH 

Capillary surfaces’ were calculated for a circular cross section with 

four rectangular reentrant notches as shown in Fig. 15A. The calculation 

was performed on 1/8 of the cross section, which is shown in Fig. 15B. 

The circle has radius 1. The notch extends inward to a depth 0.2 on the 

axis and has half-width 0.2. The corner, point G of Fig. 15B, corresponds 

to a Y of 50.77°. 

' crxt 

A series of exploratory calculations were made with 9x8 meshes 
for K=1.0 and Y=0°, 15°, 30°, 45°, and 60°. Figures 16 and 17 show 

the 9x8 meshes used for Y“60° and Y=0°* ^he case K = 1.0 and Y = 0° 
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was calculated also with the 13X12 mesh shown in Fig- 18. The cases 
with Y ^crit calculated by deleting the corner as described 

previously. 

The asymptotic solution of Eqs. (1) and (2) is not known for a corner 
which consists of a straight line intersecting a circular arc, as in Fig. 
15C. However, the arc GHE differs only slightly from its chord GJE. 

The asymptotic solution for the corner which has 'CGJE for its cylinder 
wall is the function v(p,6) defined previously. For our calculations, 
the dividing curve CDE was chosen to be a level curve of this v(p,9). 

Figures 19 and 20 show level curves of the numerical solutions for 
the 60° and 0° cases. In the 60° case, the fluid slopes gently up toward 
the cylinder wall and the corner G. In the 0° case, the fluid is drawn 
strongly up into corner G. 

Figure 21 shows the variation of the height of the capillary surface 
with the contact angle Y> various points on the cross section. The 
height is taken to be zero at the center of the circle, point 0. The 
13x12 and 9x8 meshes give somewhat different values for the surface 
height for Y^O"* The heights at points A, B, and F differ by the 
fractions 0.05/0.76, 0.02/1.52, and 0.07/1.42. The meshes do not 

give high accuracy because only a few points are used to approximate the 
corner region, where the gradient is steep. However, these calculations 
show the general nature of the solution. Greater accuracy can be obtained 
by refining the mesh and moving the level curve CDE nearer the corner G. 

Figure 22 shows the height of the surface for Y = 60° as a function 
of the arc length along the cylinder wall ABCGEF. It shows also the 
heights of the surfaces for Y = 0° and 30° along the wall from point A 
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to C and from E to F. The apparent discontinuity in the slope of the 
y=0“ surface at point B results not from an actual discontinuity, 
but because the cylinder wall changes direction by 90° at B. The 
slopes on one side and the other of the point B in Fig. 22 refer to 
two different directions. For the 60° case, the fluid rises only slightly 
from 0.305 at F to 0.506 at corner G. In the 0° and 30° cases, it rises 
steeply as G is approached. 

CDE is a level curve of v. It is a circular arc that intersects the 
cylinder wall at an angle y. The curve CDE was approximated by 4 straight 
line segments of the 9x8 mesh and by 8 segments of the 13x12 mesh. 

Table 6 shows the calculated height along these level curves. 


TABLE 6. Calculated height along CDE. 


Y 

Mesh 

u(C) 

u(D) 

u(E) 

0° 

13 X 12 

8.24 

9.83 

8.25 

15° 

9x8 

5.36 

5.76 

5.16 

30° 

9x8 

2.59 

2.65 

2.48 

45° 

9x8 

0.930 

0.933 

0.913 


The height is- quite level along CDE for the 45° case. However, 

it increases by 16% at the middle of the curve for the 0° case. For this 

case, the level curve of u is less bowed toward the corner than the 
level curve of v. The level curves of u asymptotically approach 
circular arcs as their distance to the corner approaches zero, because 
u^/v approaches 0. However, the point C is not very close to the corner. 
It is almost halfway from the corner G to point B. If a finer mesh is 

used, CDE can be chosen nearer the corner, where the level curve of v 

more closely approximates that of u. 
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Pig. 3. Height u of the capillary surface as a function of the contact 
angle y oii the edge x = 1.0 for K=0 and y=0.0, 0.6, 0.8, 
and 1.0. (XBL 775-954) 
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Pig. 4. One-fourth ellipse with a =1.0 and b = 0.6 showing values of 
(x,y) at the corners and values of L or M on. the edges. 

(XBL 775-948) 
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Fig. 10. Height u of the capillary surface at the point B as a 
function of the contact angle y for K=0 and b = 0,60, 

0.33, and 0.20 and for K=1 and b = 0.20. (XBL 776-1118) 
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Fig. 13. Height u of the capillary surface along the edge and the 

diagonal, respectively, as a function of the distance p from 
the corner. (XBL 776-1120) 
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Fig. 14. Height u of the capillary surface at the midpoint B of 

the edge of the full square as a function of contact angle Y • 

(XBL 776-1115) 
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i/8 Circle with notch 
Mesh plot y = 60° 
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Fig. 16. 9x8 mesh for a cylinder with a reentrant notch for y= 60°. 

(XBL 776-1123) 
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Fig. 17. 9x8 mesh for a cylinder with a reentrant notch for y=0°. 

(XBL 776-1122) 
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1/8 Circle with notch 
Equipotentia 1 plot K=1.0 r=0° 



Fig. 20. Level curves of the capillary surface for Y=0°* 

(XBL 776-1121) 
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Fig. 21. Height u of the capillary surface at points A, B, E, and 
F as a function of contact angle y. (XBL 775-952) 
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ABSTRACT 

In this report we show there will exist a critical contact angle 
for an elliptical cross section if the ratio b/a of the minor and 
major semiaxes is less than 0.6116. (There is no solution of the 
capillary surface Eqs. (1) and (2) for contact angles less than the 
critical angle.) We also calculate a lower bound on various 

values of b/a. These bounds appear to be close to the values of 
found by numerical solution of Eqs. (1) and (2). 
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1. DEFINITION OF THE PROBLEM 

We consider the equilihrium free surface of a liquid partly filling 
a vertical cylinder for the case in which the surface height u(x,y) is 
a single-valued smooth function of x and y, and in which there is 
sufficient liquid to cover the base of the cylinder entirely. The 
gravitational field is taken to be positive when directed vertically 
downward. The height then satisfies the equation 

V * Vu) = Ku + 2H (1) 

W = (1 + |Vu|^)^^^, 

where V is the two-dimensional operator (3/9x, 3/3y) , K = Pg/cr is the 
capillary constant, p is the difference in densities between the 
liquid and gas phases, g is the acceleration due to gravity, and O' 
is the gas-liquid surface tension. The constant 2H is determined by 
the cross-sectional shape of the cylinder, the volme of liquid, and the 
boundary condition satisfied by the free surface of the liquid at the 
cylinder wall. 

The boundary condition for a free surface that makes a contact angle 
Y with the cylinder wall is 

= cos Y at the wall, (2) 

W dn 

where 9u/3n denotes the derivative of u with respect to the outward 
directed normal at the wall. Only the case of wetting liquids 
0® < Y ^ will be considered here. (The nonwetting case can be 

obtained from it directly by means of a simple transformation.) 
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2. CRITERION FOR THE EXISTENCE OF A CRITICAL CONTACT ANGLE FOR 
ELLIPTICAL CROSS SECTIONS 

Let A denote the area and L the perimeter of the cross section 
of a cylinder. The following theorem holds [1,2]. 

Theorem . Suppose there is a point p on the boundary at which the 
curvature is greater than L/A; then there exists a critical contact 
angle such that there is no solution of Eqs . (1) and (2) in a neighbor- 
hood of p for 0 < Y Y 

^ crit 

We shall apply this theorem to an elliptical cross section. For 
this application we m\ist calculate the area and perimeter of an ellipse 
and the curvature at any point on it. The area- of an ellipse is 

A = irab (3) 

where a and b are the semimajor and semiminor axes. 

The ellipse is described by the equation 

y(x) = b Ci - x^/a^)^^^ . (4) 


The perimeter of the ellipse is 

a 


L = 4 J* [1 + dx 


a n 1 , 2 . 2 , 2 , 24/2 

- 4 r [1- Cl-b M ^ /a ] ■ dx 

~ ^ n 2 , 2 . 1/2 *. 

o (1 - X /a ) 


Let x/a. = sin 9, and 


1 1,2/ 2 

m = 1 ~ b /a 


(5) 


Then 


L 


4aE (m) 




( 6 ) 
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where 


E(m) = 




(1 - m sin^ 9)^ d0 


(7) 


is the complete elliptic integral of the second kind. Values of L/A 
for a = 1 and various values of b are given in Table I to the 
indicated number of decimal places . 


Table I. Values of L/A for a = 1 and various values of b. 


b 

m 

L/A 

0.2 

0.96 

6.688 

0.25 

0.9375 

5.461 

0.33 

0.8911 

4.290 

0.4 

0.84 

3.663 

0.5 

0.75 

3.084 

0.6 

0.64 

2.708 

These will be used in a 

later section of 

this report. 


Next we shall calculate the curvature. The differential of arc 
length is 

ds = [1 + dx 

^ A 

The unxt vector t tangent to the curve y(x) in two dimensions 


is 


t = (dx/ds, dy/ds) 

= [1 + 


(i.yO . 


( 8 ) 
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The curvature C and the mit vector fi normal to the curve y(x) 
are defined by 

Cn = dt/ds 

dt/ds = [1 + (y*)^] d£/dx. 

= y" II + (y’)^r^ ,(-yM) . 


(9) 


Therefore 


and 


n = sgn (y”) [1 + (y’)^3 (-y’,1) 

C = |y'* [1 + (y')^r^^^l . 


( 10 ) 

( 11 ) 


Note that C also can be written 

C 

Let y(x) be an ellipse. Then 


■\^i 

' dx ) 


y' [1 + (y')^3^^^ 




( 12 ) 


Thus 


y* = -b^x/a^y 

II 1.^/ 2 3 
yi I = _b /ay 


. 4 -4,, 4 2 , . 4 2.3/2 

C = a b /(a y + b x ) 


The curvature is maximum at the point x = a, y==0 where the major 
axis intersects the ellipse. Thus 

2 


C = a/b 
max 


(13) 


According to the theorem, there will exist a critical contact angle for 

ellipses with C greater than L/A, that is, for 
max 

b/a < 7r/4E(m) . (14) 

This inequality can be solved numerically. The result is that a 
critical angle will exist for ellipses with b/a < 0.6116. 
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3. A LOWER BOUND ON Y • 

crxt 

Consider the cross section of the cylinder shown in Fig. 1. Let 
A be the area of the entire cross section and L be its perimeter. 

The cross section is cut by a curve F, which Intersects the perimeter 
at points and ^ 2 ^ whose length we also denote by F. Let 

A* denote the part of the domain cut off by F and denote also the 
area of this part. Let L* denote the part of the perimeter cut off by 
F and denote also the length of this part. 

If a solution to Eqs . (1) and (2) exists for K = 0 and contact 
angle y, then integration of Eq. (1) over A* gives 

2HA* = / V • (“ Vu) dA = / d Vu)»n dL 

A* ^ F+L* ^ 

= L* cos y + S (.h Vu)*n dL . 

■ F ” 

Now |h*Vu| therefore 

-F < / (~ Vu) • n dL < r . 

F ^ 

Thus if a solution exists for contact angle y, then for any 
curve F 

-r < 2HA* - L* cos y < F . (15) 


This can be written 

2HA* - r ^ ^ 2HA* + r 

T-T cos y ^ r-r . 

L* L* 

Integration of Eq. (1) over the entire cross section gives 


( 16 ) 


2HA = L cos Y . 
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Eliminating H from Eq, (15) gives 

-r •< (L* - A*L/A) cos Y ^ r • 

For 0 ^ Y < 90° and L* - A*L/A nonzero, this can be written 

cos Y j 

where 

^ " |L* - A*L/A| ■ 

Equations (16) and (17) are Lemma 6 and Corollary 3.2 of Eef. 2 for the 
special case of a two-dimensional domain. 

Equation (17) holds for each contact angle for which a solution to 
Eqs. (1) and (2) exists. In particular it holds for Let 

be the minimum of V with respect to all possible curves F. is 

an upper bound for cos 


cos y . V . 

crxt o 

This gives a lower bound on Y • • 

^ ' crxt 


(19) 


4. MINIMI2ATI0N OF V FOR FIXED p^ AND 

The minixnization of V can be done in two steps. First, find F 
that minimizes V for a fixed pair of intersection points p^ and P 2 S 
then vary p^ and p^. It can be shown that the minimizing F is the 
arc of a circle [2]. 

Let R be the radius of the circular arc F, and let 0 be the 
center of the circle as shown in Fig. 2. Let 2y be the length of the 
chord from p^ to P 2 * Let A be the area of the region bounded by 
this chord and by F. Let 20 be the angle 0 p 2 - Then the follow- 
ing, relations hold. 



0 
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( 20 ) 

( 21 ) 

( 22 ) 


= arcsin (y/R) 

r = 2R0 

2 

= R 0 - yR cos 6 

Let us hold and P 2 fixed and vary R. Then y is fixed* but 6, 

r, and vary. Also 5 a* = - 5A^ because A* + A^^ is fixed-. 

Varying R in Eqs . (20)-(22) and eliminating 60 gives 

6r = 2(0 - tan 6) 6R 

6A* = -2R(6 - tan 0) 6R 

thiis 

6A* = -R 6r . . (23) 

Let (L* - A*L/A) be positive; then 

(L* - A*l/A) 6V = 6r + (VL/A) 6A* . (24) 

Let Vp be the minimum value of V for fixed points p^^ and P 2 > 
and let R^ be the minimizing radius. Setting 6V = 0 in Eq. (24) and 
using Eq. (23) gives 

.R^ = A/(LV ) . (25) 

P P 

5. MINIMIZATION OF V WITH AND P 2 

Next, let us hold R fixed and vary p^^ and P 2 * Let <j>^ denote 
the acute angle of Intersection of T with the cross section at 
and let ^>2 <ieJ^ote the corresponding angle at P 2 * Let 6 L*j^ be the 
variation in L* at p^, and ^L 2 * be the variation at P 2 * Then 

61 = cos (|)^ ^^ 1 * ^2 ^^2* 

6A* = (r/2) sin <J)^ 5L^* + (F/2) sin <SL 2 * 


(27) 
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Let (L* - A*L/A) be positive; then 

(L* - A*L/A) 6v = 6r + (YL/A) 6 A* - V 6L* (28) 

Setting 5V = 0 in Eq. (28) and using Eqs. (26) and (27) gives two 
equations for the coefficients of the independent variations 6L^* 
and 

cos + U sin ^ (29) 

cos 4>2 U sin ({>2 = V , (30) 

where 

TJ = VLr/2A 

For a given value of V, L, F, and A, Eq. (29) has the solutions 

sin (|)^ = |u + + (1 + U^)(l - V^)]^^^j“/(1 + U^) . 

For V < 1, Eq. (29) has one positive solution for sin (j>^. For 

V > 1, it has two positive solutions; however, V > 1 gives no useful 

bound on cos Y ■ , so we shall ignore this case. Thus if V has a 
'crxt ^ 

relative minimuiii that is less than one, then mil equal (j >2 for 

the minimizing T. 

Consider a cross section that is symmetric about some straight 
line, as is the ellipse of Fig. 3 about the x axis. Let the curve 
r have intersection points p^ and P 2 on opposite sides of that line. 
If V has a relative minimum that is less than one, then the minimizing 
p^ and P 2 will be symmetric about that line. 
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6. LOWER BOUND ON Y • 5'OR ELLIPTICAL CROSS SECTIONS 

crit 

We apply the preceding results to the elliptical cross section 
shown in Fig. 3. The ellipse is 


/ ^ 2 , 2 . 1/2 

y(x) = b(l“x/a) 


(31) 


Let p^ be the point [x^, y(x^)]; it is sufficient to take to be 

[x^, - y(x^)]. Let be the area + A*; then 


A 


2 


2/ y(x) 

X, 


dx 


(32) 


A* = A^ “ A^ (33) 

^ 9 1 /9 

L* = 2 j* [1 + (y')^]^^^ dx . (34) 


Equations (18), (20)-(22), (25), and (31)-(34) can be solved for 

for each value of x^^. By calculating for a sequence of values of 

X- , we find the minimum value V . 

1* o 

The circular arc T must lie inside the ellipse. This is assured 
if R is greater than R^, R^ is the length of the line 0^ that 
is perpendicular to the ellipse at p^ as shown in Fig. 3. 


^1 = 

[y(xj_) 3^ + b^x^ / a^ 

(35) 

R = 

max (R, , R ) 

1’ p 

(36) 


The preceding equations give the following results for the lower 

bound on Y - a = 1 and various values of b . 

' crtt 
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b 

letter bound 

0.20 

35.12“ 

0.25 

26.95° 

0.33 

16.88° 

0.40 

10.32° 

0.50 

3.71° 

0.60 

0.12° 


These bounds appear to be close to the values of found by 
numerical solution of Eqs. (1) and (2) 13]. The lower bound as a 
function of b is shown in Fig. 4. 
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Fig. 1. Cross section of a cylinder. (XBL 776-949) 



Fig. 2. Cross section of a cylinder with variables defined. (XBL 776-949) 



Fig. 3 


Cross section of an ellipse. (XBL 776-949) 
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Fig. 4 



. Lower bound on Y . as a function of th.e ratio b of the 

crit 

semiminor and semimajor axes. (XBL 776-1117) 
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